Comparison  of  Optimal  Design  Methods  in  Inverse  Problems 


H.  T.  Banks1,2,  Kathleen  Holm1,2  and  Franz  Kappel2,3 


1  Center  for  Quantitative  Sciences  in  Biomedicine 
North  Carolina  State  University 
Raleigh,  NC  27695-8213 

2Center  for  Research  in  Scientific  Computation 
North  Carolina  State  University 
Raleigh,  NC  27695-8212 
and 

institute  for  Mathematics  and  Scientific  Computation 
University  of  Graz 
Graz,  Austria  A8010 

May  11,  2011 


Abstract 

Typical  optimal  design  methods  for  inverse  or  parameter  estimation  problems  are  designed 
to  choose  optimal  sampling  distributions  through  minimization  of  a  specific  cost  function  related 
to  the  resulting  error  in  parameter  estimates.  It  is  hoped  that  the  inverse  problem  will  produce 
parameter  estimates  with  increased  accuracy  using  data  from  the  optimal  sampling  distribution. 
We  present  a  new  Prohorov  metric  based  theoretical  framework  that  permits  one  to  treat  suc¬ 
cinctly  and  rigorously  any  optimal  design  criteria  based  on  the  Fisher  Information  Matrix  (FIM). 
A  fundamental  approximation  theory  is  also  included  in  this  framework.  A  new  optimal  de¬ 
sign,  S'.E-optimal  design  ( standard  error  optimal  design),  is  then  introduced  in  the  context  of 
this  framework.  We  compare  this  new  design  criteria  with  the  more  traditional  D-optimal  and 
£,-optinral  designs.  The  optimal  sampling  distributions  from  each  design  are  used  to  compute 
and  compare  standard  errors;  the  standard  errors  for  parameters  are  computed  using  asymptotic 
theory  or  bootstrapping  and  the  optimal  mesh.  We  use  three  examples  to  illustrate  ideas:  the 
Verhulst-Pearl  logistic  population  model  [7],  the  standard  harmonic  oscillator  model  [7]  and  a 
popular  glucose  regulation  model  [10,  13,  21]. 
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1  Introduction 


Mathematical  models  are  used  to  describe  dynamics  arising  from  biological,  physical  and  engi¬ 
neering  systems.  If  the  parameters  in  the  model  are  known,  the  model  can  be  used  for  simulation, 
prediction,  control  design,  etc.  However,  typically  one  does  not  have  accurate  values  for  the  pa¬ 
rameters.  Instead,  one  must  estimate  the  parameters  using  experimental  data.  The  simulation  and 
predictive  capabilities  of  the  model  depend  on  the  accuracy  of  the  parameter  estimates.  A  major 
question  that  experimentalists  and  inverse  problem  investigators  alike  often  face  is  how  to  best  col¬ 
lect  the  data  to  enable  one  to  efficiently  and  accurately  estimate  model  parameters.  This  is  the 
well-known  and  widely  studied  optimal  design  problem. 

Traditional  optimal  design  methods  (D-optimal,  E-optimal,  c-optinral)  [1,  8,  14,  15]  use  infor¬ 
mation  from  the  model  to  find  the  sampling  distribution  or  mesh  for  the  observation  times  (and/or 
locations  in  spatially  distributed  problems)  that  minimizes  a  design  criterion,  quite  often  a  function 
of  the  Fisher  Information  Matrix  (FIM).  Experimental  data  taken  on  this  optimal  mesh  is  then 
expected  to  result  in  accurate  parameter  estimates. 

Here  we  formulate  the  classical  optimal  design  problem  in  the  context  of  general  optimization 
problems  over  distributions  of  sampling  times.  We  present  a  new  Prohorov  metric  based  theoretical 
framework  that  allows  one  to  treat  succinctly  and  rigorously  any  optimal  design  criteria  based  on 
the  FIM.  A  fundamental  approximation  theory  is  also  included  in  this  framework.  A  new  optimal 
design,  S' D-optimal  design  ( standard  error  optimal  design),  is  then  introduced  in  the  context  of 
this  framework.  We  compare  this  new  design  criteria  with  the  more  traditional  D-optimal  and 
D-optimal  designs.  We  consider  the  performance  of  these  three  different  optimal  design  methods 
for  three  different  dynamical  systems:  the  Verhulst-Pearl  logistic  population  model,  a  harmonic 
oscillator  model  and  a  simple  glucose  regulation  model.  DD-optimal  design  was  first  introduced  in 
[5].  The  goal  of  SD-optimal  design  is  to  find  the  observation  times  r  =  {t,}  that  minimize  the 
sum  of  squared  normalized  standard  errors  of  the  estimated  parameters  as  defined  by  asymptotic 
distribution  results  from  statistical  theories  [4,  6,  12,  20].  D-optimal  and  D-optimal  design  methods 
minimize  functions  of  the  covariance  in  the  parameter  estimates  [1,  8,  15].  D-optimal  design  finds 
the  mesh  that  minimizes  the  volume  of  the  confidence  interval  ellipsoid  of  the  asymptotic  covariance 
matrix.  D-optimal  design  minimizes  the  largest  principle  axis  of  the  confidence  interval  ellipsoid  of 
the  asymptotic  covariance  matrix. 

In  an  effort  to  provide  a  reasonably  fair  comparison,  for  each  optimal  design  method,  standard 
errors  are  computed  by  several  methods  using  the  optimal  mesh.  The  optimal  design  methods  are 
compared  based  on  these  standard  errors.  Not  surprisingly,  we  find  that  5D-optimal  design  often 
results  in  smaller  standard  errors  compared  with  the  other  optimal  design  method;  this  is  likely  be¬ 
cause  S’D-optimal  design  optimizes  directly  on  the  standard  errors  themselves  while  the  D-optimal 
and  D-optimal  methods  minimize  other  functions  related  to  the  standard  errors  through  the  FIM. 


2  Optimal  Design  Formulations 

Following  [5],  we  introduce  a  formulation  of  ideal  inverse  problems  in  which  continuous  in  time 
observations  are  available-while  not  practical,  the  associated  considerations  provide  valuable  insight. 
A  major  question  in  this  context  is  how  to  choose  sampling  distributions  in  an  intelligent  manner. 
Indeed,  this  is  the  fundamental  question  treated  in  the  optimal  design  literature  and  methodology. 
Underlying  our  considerations  is  a  mathematical  model 
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x{t)  =g(t,x(t),q), 
x(0)  =x0, 

f(t,6)  =  C(x(t,0)),  iG[0,T] 


(1) 


where  x[t )  G  is  the  vector  of  state  variables  of  the  system,  f(t,  8)  G  Mm  is  the  vector  of  observable 
or  measurable  outputs,  q  G  Mr  are  the  system  parameters,  8  =  (q,x o)  G  Wp,p  =  r  +  n  is  the  vector 
of  system  parameters  plus  initial  conditions  xq,  while  g  and  C  are  mappings  jj1+ri+r  — Mn  and 
Mn  — >■  Rm,  respectively.  To  consider  measures  of  uncertainty  in  estimated  parameters  [4],  one  also 
requires  a  statistical  model.  Our  statistical  model  is  given  by  the  stochastic  process 

Y(t)  =  f(t,e0)+S(t).  (2) 

Here  £  is  a  noisy  random  process  representing  measurement  errors  and,  as  usual  in  statistical  for¬ 
mulations  [4,  5,  20],  8q  is  a  hypothesized  “true”  value  of  the  unknown  parameters.  We  make  the 
following  standard  assumptions  on  the  random  variable  £(t): 

E(£(£))  =  0,  i  G  [0,T], 

Var  £(t)  =  ctq,  £g[0,T], 

Co v(£(t)£(s))  =  (Tq 5(t  —  s),  t,s  G  [0,T], 

where  5(s)  =  1  for  s  =  0  and  5(s )  =  0  for  s  /  0.  A  realization  of  the  observation  process  is  given  by 

y(t)  =  f(t,0o)  +  e(t)y  t  G  [0,  T], 
where  the  measurement  error  e(t)  is  a  realization  of  £{t). 

We  introduce  a  generalized  weighted  least  squares  criterion 

J(y^9)=[  — 77T2  {v(t)  ~  f(t.,6))2  dP(t),  (3) 

Jo  a\t) 

where  P  is  a  general  measure  on  [0,T].  We  seek  the  parameter  estimate  8  by  minimizing  J(y,8 )  for 
8.  Since  P  represents  a  weighting  of  the  difference  between  data  and  model  output,  we  can,  without 
loss  of  generality,  assume  that  P  is  a  bounded  measure  on  [0,T], 

If,  for  points  r  =  {tj},  t\  <  ■  ■  ■  <  ijv  in  [0,  T],  we  take 

N 

Pr=YlAti,  (4) 

2=1 

where  Aa  denotes  the  Dirac  delta  distribution  with  atom  {a},  we  obtain 

N  1 

Jd{y,8)  =  ^j——-1{y{ti)  -  f{th8)f ,  (5) 

*= 1 

which  is  the  weighted  least  squares  cost  functional  for  the  case  where  we  take  a  finite  number  of 
measurements  in  [0,  T].  Of  course,  the  introduction  of  the  measure  P  allows  us  to  change  the  weights 
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in  (5)  or  the  weighting  function  in  (3).  For  instance,  if  P  is  absolutely  continuous  with  density  m(-) 
the  error  functional  (3)  is  just  the  weighted  L2-norm  of  y(-)  —  with  weight  m(-)/cr(-)2. 

To  facilitate  our  discussions  we  introduce  the  Generalized  Fisher  Information  Matrix  (GFIM) 

F(P,e  0)=  [T^—\/Jf(s,e0)Vef(s,e0)dP(s),  (6) 

Jo  (s) 

where  V#  is  a  row  vector  given  by  (d$1, . . .  ,dgp )  and  hence  Vgf  is  an  m  x  p  matrix.  It  follows  that 
the  usual  discrete  FIM  corresponding  to  PT  as  in  (4)  is  given  by 

N  1 

F(t)  =  F(Pr,  do)  =  ]T  ——Wefitj,  e0 )TV0/(tJ,  do).  (7) 

j=1  °  V-j) 

Subsequently  we  simplify  notation  and  use  r  =  {t-i}  to  represent  the  dependence  on  P  =  PT  when  it 
has  the  form  (4).  When  one  chooses  P  as  simple  Lebesgue  measure  then  the  GFIM  reduces  to  the 
continuous  FIM 


Fc=  [T  -4r^vof(s,e0)TVef(s,d0)ds.  (8) 

Jo  ^(s) 

The  major  question  in  optimal  design  of  experiments  is  how  to  best  choose  P  in  some  family  P(0,T) 
of  observation  distributions.  We  observe  that  one  optimal  design  formulation  we  might  employ  is  a 
criterion  that  chooses  the  times  r  =  {t;}  for  PT  in  (6)  so  that  (7)  best  approximates  (8)-i.e.,  one 
minimizes  | Fq  —  T(r)|  over  r  where  |  •  |  is  the  norm  in  Wpxp-see  [5].  We  do  not  consider  this  design 
here,  but  rather  focus  on  the  SE-optimal  design  also  proposed  in  [5]  and  its  comparison  to  more 
traditional  designs. 

The  introduction  of  the  measure  P  above  allows  for  a  unified  framework  for  optimal  design 
criteria  which  incorporates  all  the  popular  design  criteria  mentioned  in  the  introduction.  As  already 
noted,  the  GFIM  F(P,9 )  introduced  in  (6)  depends  critically  on  the  measure  P.  We  also  remark 
that  we  can,  without  loss  of  generality,  further  restrict  ourselves  to  probability  measures  on  [0, T\. 
Thus,  let  P(0,T)  denote  the  set  of  all  probability  measures  on  [0,  T]  and  assume  that  a  functional 
J  :  M.pxp  — >  M+  of  the  GFIM  is  given.  The  optimal  design  problem  associated  with  J  is  one  of 
finding  a  probability  measure  P  €  V(0,T)  such  that 

J{F(P,6  0))=  min  j(F(P,0o)).  (9) 

A  general  theoretical  framework  for  existence  and  approximation  in  the  context  of  V(0 ,T)  taken 
with  the  Prohorov  metric  [2,  11,  16,  19]  is  given  for  these  problems  in  Section  4  of  [5].  In  particular, 
this  theory  permits  development  of  computational  methods  using  weighted  discrete  measures  (i.e. , 
weighted  versions  of  (4)). 

2.1  Theoretical  Summary 

To  summarize  and  further  develop  the  theoretical  considerations  that  are  the  basis  of  our  efforts 
here,  we  first  recall  that  the  Prohorov  metric  p  on  the  space  V(0,T)  of  probability  measures  on 
the  Borel  subsets  of  [0,T]  can  be  defined  [2,  11,  16,  19]  in  terms  of  probabilities  on  closed  subsets 
of  [0,  T]  and  their  neighborhoods.  However  for  our  uses  here  it  is  far  more  useful  to  work  with  an 
equivalent  characterization  in  terms  of  convergences  when  viewing  the  probability  measures  ,  T) 
as  a  subset  of  the  topological  dual  Cb[ 0,  T]*  of  the  bounded  continuous  functions  on  [0,  T]  taken  with 
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the  supremum  norm.  More  precisely,  p-convergence  is  equivalent  to  weak*-convergence  on  V(0 ,T) 
when  considering  V(0,T)  as  a  subset  of  Cb[0,T]*.  It  is  then  known  that  (P(0,  T),p)  is  a  complete, 
compact  and  separable  metric  space.  (We  will  hereafter  just  denote  this  space  by  P(0,  T)  since  the 
p  will  be  understood.) 

Our  first  observation  is  that  the  GFIM  as  defined  in  (6)  is  p  continuous  on  P(0,  T )  for  problems 
in  which  the  observation  functions  /(•,#)  are  continuously  differentiable  on  [0,T],  Thus,  whenever 
J  :  Mpxp  — >  R+  is  continuous  we  find  that  P  — >  J(F(P ,  6))  is  continuous  from  P(0,  T)  to  R+.  Since 
P(0,  T )  is  p  compact,  we  obtain  immediately  the  existence  of  solutions  for  the  optimization  problems 

Pj  =  argminpep(0T)l7  (F(P,6q)).  (10) 

Our  second  observation  is  related  to  the  separability  of  P(0,  T)  and  in  particular  to  the  den¬ 
sity  of  finite  convex  combinations  over  rational  coefficients  of  Dirac  measures  Aa  with  atoms  at  a. 
Specifically,  one  can  prove  [2]  that  the  set 

k  k 

V0(0,T)  :=  jp  6  P(0,T)  P  =  '^2pjAtj,  k  <G  N+,  tj  <G  7o,  pj  >  0,  pj  rational,  ^ pj  =  lj 

j= i  i= 1 

is  dense  in  P(0 ,T)  in  the  Prohorov  metric  p.  Here  7o  =  }^=i  is  a  countable,  dense  subset  of 

[0,  T].  In  short,  the  set  of  P  €  P(0,  T)  with  finite  support  in  7o  and  rational  masses  is  dense  in 
P(0 ,T).  This  leads,  for  a  given  choice  J ,  to  approximation  schemes  for  Pj  as  defined  in  (10).  To 
implement  these  for  a  given  choice  of  J  (examples  are  discussed  below)  would  require  approximation 
by  Pjp  t  y  =  Pj  A-tj  in  the  GFIM  (6)  and  then  optimization  over  appropriate  sets  of  {pj ,  t3 } 

in  (10)  with  P  replaced  by  P^  t  y.  For  a  fixed  N.  existence  of  minima  in  these  problems  follow 
from  the  theory  outlined  above.  In  standard  optimal  designs  these  problems  are  approximated  even 
further  by  fixing  the  weights  or  masses  pj  as  pj  =  -j^  (which  then  becomes  simply  a  scale  factor  in  the 
sum)  and  searching  over  the  {tj}.  This,  of  course,  is  equivalent  to  replacing  the  Pj^  t  ,  by  PT  of  (4) 
in  (6)  and  searching  over  the  r  =  {tj}  for  a  fixed  number  N  of  grid  points.  This  embodies  the  tacit 
assumption  of  equal  value  of  the  observations  at  each  of  the  times  {tj}.  We  observe  that  weighting 
of  information  at  each  of  the  observation  times  is  carried  out  in  the  inverse  problems  via  the  weights 
( j(tj )  for  observation  variances  in  (5).  We  further  observe  that  the  weights  {pj}  in  Pj^  t  ,  are  related 
to  the  value  of  the  observations  as  a  function  of  the  model  sensitivities  Vgf(tj,  do)  in  the  FIM  while 
the  weights  are  related  to  the  reliability  in  the  data  measurement  processes.  We  note  that  all 

of  our  remarks  on  theory  related  to  existence  above  in  the  general  probability  measure  case  also  hold 
for  this  discrete  minimization  case. 

The  formulation  (10)  incorporates  all  strategies  for  optimal  design  which  entail  optimization 
of  a  functional  depending  continuously  on  the  elements  of  the  Fisher  information  matrix.  In  case 
of  the  traditional  design  criteria  mentioned  in  the  introduction,  J  is  the  determinant  (D-optimal), 
the  smallest  eigenvalue  (E-optimal),  or  a  quadratic  form  (c-optimal),  respectively,  of  the  inverse  of 
the  Fisher  information  matrix.  Specifically,  the  optimal  design  methods  we  consider  are  S'P-optimal 
design,  D-optimal  design,  and  D-optimal  design.  The  design  cost  functional  for  the  DD-optimal 
design  method  is  given  by  (see  [5]) 

=izw(FP“’ 

i=l  °0,i 

where  F  =  P(r)  is  the  FIM,  defined  above  in  (7),  9 o  is  the  true  parameter  vector,  and  p  is  the 
number  of  parameters  to  be  estimated.  Note  that  both  inversion  and  taking  the  trace  of  a  matrix  are 
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continuous  operations.  We  observe  that  T)”1  =  SE,(0q)2.  Therefore,  iSlE-optimal  design  minimizes 
the  sum  of  squared  normalized  standard  errors. 

id-optimal  design  minimizes  the  volume  of  the  confidence  interval  ellipsoid  for  the  covariance 
matrix  (£q  =  F~ 1).  The  design  cost  functional  for  id-optimal  design  is  given  by  (see  [8,  15]) 

J(F)  =  det  (F-1). 

Again  we  note  that  taking  the  determinant  is  a  continuous  operation  on  matrices  so  that  Jr>  is 
continuous  in  F  as  required  by  the  theory. 

E-optimal  design  minimizes  the  principle  axis  of  the  confidence  interval  ellipsoid  of  the  covariance 
matrix  (defined  in  the  asymptotic  theory  summarized  in  the  next  section).  The  design  cost  functional 
for  id-optimal  design  is  given  by  (see  [1,  8]) 

J(F)  =  max—, 

^ i 

where  An  i  =  1 . . . p  are  the  eigenvalues  of  F.  Therefore  4~,  i  =  1 ..  .p,  corresponds  to  the  eigenvalues 
of  the  asymptotic  covariance  matrix  =  F~l. 


2.2  Constrained  Optimization  and  Implementations 

Each  optimal  design  computational  method  we  employ  is  based  on  constrained  optimization  to 
find  the  mesh  of  time  points  r*  =  {t*},  i  =  1 , ...  ,1V  that  satisfy 

J(F(T*,6o,))  =  mmJ(F(T,0o)), 
tE  / 

where  T  is  the  set  of  all  time  meshes  such  that  0  <  t\  <  •  •  •  <  ti  <  tj+i  <  •  •  •  <  tjy  <  T. 

These  optimal  design  methods  were  implemented  using  constrained  optimization  algorithms, 
either  MATLAB’s  fmincon  or  SolvOpt,  developed  by  A.  Kuntsevich  and  F.  Kappel  [17],  with  four 
variations  on  the  constraint  implementation.  We  denote  these  different  constraint  implementations 
(which  result  in  different  parameter  and  SE  outcomes  even  in  cases  where  the  {t,}  are  initially 
required  to  satisfy  similar  constraints)  by  (Cl)  —  (C4).  Complete  details  of  the  differences  in  the 
algorithms  are  given  in  an  appendix. 

(Cl)  The  first  constraint  implementation  on  the  time  points  is  given  by,  t\  >  0,  tjv  <  T  and  ti  <  t?.+i, 
such  that  the  optimal  mesh  may  or  may  not  contain  0  and  T.  In  this  case  we  optimize  over  N 
variables. 

(C2)  The  second  constraint  implementation  is  carried  out  in  the  same  manner  as  the  first,  except 
that  the  optimal  mesh  contains  0  and  T.  Hence  we  effectively  optimize  over  N  —  2  variables. 

(C 3)  The  third  constraint  implementation  on  the  time  points  is  given  by  ti  =  C-i  +  ty,  i  =  2, . . . ,  N  — 
1,  t\  =  0  and  tjsi  =  T,  with  >  0,  i  =  2, . . . ,  N  —  1,  and  v-2  +  . . .  +  un~i  <  T.  Note  that 
the  optimal  mesh  always  contains  0  and  T  as  we  optimize  over  N  —  2  variables  using  slightly 
different  inequality  constraints. 

(C4)  The  last  constraint  implementation  on  the  time  points  is  given  by,  ti  =  +  zzj,  i  =  2, . . . ,  N, 

and  t\  =  0  with  >  0,  i  =  2, . . . ,  N,  and  z/2  +  . . .  +  =  T.  This  constraint  is  implemented  by 

defining  z/jv  =  T  —  ui-  r^^Le  opthnal  mesh  again  contains  0  and  T,  and  we  also  optimize 

over  N  —  2  variables  but  an  equality  constraint  is  added  to  the  constraint  system. 
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3  Standard  Error  Methodology 

We  begin  by  finding  the  optimal  discrete  sampling  distribution  of  time  points  r  =  {U}^L1,  for 
a  fixed  number  N  of  points  in  a  fixed  interval  [0,T],  using  one  of  three  optimal  design  methods 
described  above.  These  three  optimal  design  methods  are  then  compared  based  on  the  standard  er¬ 
rors  computed  for  parameters  using  these  sampling  times.  Since  there  are  different  ways  to  compute 
standard  errors,  we  will  compare  the  optimal  design  method  using  different  techniques  for  computing 
the  standard  errors.  In  the  following  sections  we  will  describe  the  methods  for  computing  standard 
errors.  First  we  consider  the  scalar  observation  case  ( m  =  1). 


3.1  Asymptotic  Theory  for  Computing  Standard  Errors 

Once  we  have  an  optimal  distribution  of  time  points  we  will  obtain  data  or  simulated  data, 
{Ui}iLii  a  realization  of  the  random  process  {Tj}^1,  corresponding  to  the  optimal  time  points, 
r  =  {ti}fL1.  Parameters  are  then  estimated  using  inverse  problem  formulations  as  described  in  [4]. 
Since  the  variance  Var(£(t))  =  Uq  is  assumed  to  be  constant,  the  inverse  problem  is  formulated  using 
ordinary  least  squares  (OLS).  The  OLS  estimator  is  defined  by 


N 

©ols  =  ©ols  =  argmin^]^  -  f(tj,9)}2. 

3= 1 

The  estimate  6*ols  is  defined  as 

N 

0OLS  =  ^ols  =  argmm  ~ 

3= 1 


To  compute  the  standard  errors  of  the  estimated  parameters,  we  first  must  compute  the  sensi¬ 
tivity  matrix 

d[Cx[tj))  _  df(tj,d ) 


Xj,k  ~ 


d9u 


d9u 


-,  for  j  =  1,  . . .  N,  k  =  1, . .  .p. 


Note  that  x  =  a  is  an  N  x  p  matrix.  The  true  constant  variance 


( 7 


=  — E 
N 


N 


o)]2 


can  be  estimated  by 

1  N 

^OLS  =  Jf—  ~  ^0LS)]J- 

1  3=1 

The  true  covariance  matrix  is  approximately  (asymptotically  as  N  — >•  oo)  given  by, 


~  °o[xT (Qo)x(8o)]  1- 


Note  that  the  approximate  Fisher  Information  Matrix  (FIM)  is  defined  by 


F  =  F(T)  =  F(r,d  o)  =  (S^)-1, 


(11) 
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and  is  explicitly  dependent  on  the  sampling  times  r. 

When  the  true  values,  9q  and  a2,  are  unknown,  the  covariance  matrix  is  estimated  by 

ps  £N0ols)  =  o-ols[xT(^ols)x(^ols)]_1-  (12) 

The  corresponding  FIM  can  be  estimated  by 

F{t)  =  Hr,  Oqls)  =  (X^ols))-1.  (13) 

The  asymptotic  standard  errors  are  given  by 

SEk(90)  =  )J(Zg)kk,  k  =  1, . . .  ,p.  (14) 

These  standard  errors  are  estimated  in  practice  (when  9q  and  <7o  are  not  known)  by 

SEk(9 ols)  =  \J 0OLs))kk,  k  =  l,...,p.  (15) 

It  can  be  shown,  under  certain  conditions,  for  N  — >  oo,  that  the  estimator  0qLS  is  asymptotically 
normal  [20];  i.e. ,  for  N  large 

©ols-A/^o,^).  (16) 

3.2  Monte  Carlo  Method  for  Asymptotic  Standard  Errors 

To  account  for  the  variability  in  the  asymptotic  standard  errors  due  to  the  variability  in  the 
residual  errors  in  the  simulated  data,  we  use  Monte  Carlo  trials  to  examine  the  average  behavior. 
For  a  single  Monte  Carlo  trial,  we  generate  simulated  data  on  the  optimal  mesh  {tj}jL i, 

Vj  =  f(tj,6o)  +  ej,  j  =  1, ...  IV, 

where  ej  are  realizations  of  £j  ~  Af( 0,  a2)  for  j  =  1, . . . ,  N.  Parameters  are  estimated  using  OLS  with 
initial  parameter  guess  9 0  =  1.4#o,  where  9q  are  the  true  parameters.  Standard  errors  are  estimated 
using  asymptotic  theory  (15).  The  parameter  estimates  and  their  estimated  standard  errors  are 
stored,  and  the  process  is  repeated  with  new  simulated  data  corresponding  to  the  optimal  mesh  for 
M  =  1000  Monte  Carlo  trials.  The  average  of  the  M  =  1000  parameter  estimates  and  standard 
errors  are  used  to  compare  the  optimal  design  methods  in  one  of  our  examples. 

3.3  The  Bootstrapping  Method 

An  alternative  way  of  computing  parameter  estimates  and  standard  errors  uses  the  bootstrapping 
method  [6].  Again  we  outline  this  for  the  case  of  scalar  ( m  =  1)  observations. 

As  in  the  previous  section,  assume  we  are  given  experimental  data  (r/i,  t\), . . . ,  ( yjv,  Cv )  from  the 
following  underlying  observation  process 


Yj  =  f(tj,90)  +  £j,  (17) 

where  j  =  1  and  the  £j  are  independent  identically  distributed  (iid)  from  a  distribution 

E  with  mean  zero  (E (£j)  =  0)  and  constant  variance  <7g,  and  9q  is  the  “true”  parameter  value. 
Associated  corresponding  realizations  of  Yj  are  given  by 


Vj  =  f(tj,90 )  +  ej. 


The  bootstrapping  algorithm  is  presented  for  sample  points  corresponding  to  the  tj,  j  =  1 ...  N. 
To  compare  the  optimal  design  methods  based  on  their  bootstrapping  standard  errors,  we  will  take 
our  sample  points  corresponding  to  the  optimal  time  distribution  (r  =  {U}^).  The  different  optimal 
design  methods  are  described  below. 

The  following  algorithm  [6]  can  be  used  to  compute  the  bootstrapping  estimate  Oboot  of  Oq  and 
its  empirical  distribution. 

1.  First  estimate  9°  from  the  entire  sample,  using  OLS. 

2.  Using  this  estimate  define  the  standardized  residuals: 

fs  =  \lwh) 

for  j  =  1  ,...,N.  Then  {fi,. . .  ,vn}  are  realizations  of  iid  random  variables  Rj  from  the 
empirical  distribution  Fn,  and  p  for  the  number  of  parameters.  Observe  that 

N  N 

E(fj  I-Fjv)  =  N-1  J2  O  =  0,  Var(fj \FN)  =  JV"1  f2  =  d2. 

3= 1  3= 1 

Set  m  =  0. 

3.  Create  a  bootstrap  sample  of  size  N  using  random  sampling  with  replacement  from  the  data 
(realizations)  {fi,. . .  Fn}  to  form  a  bootstrap  sample  {r™, . . . ,  r\ y}. 

4.  Create  bootstrap  sample  points 

y™  =  ntj,e°)  +  r™, 

where  j  =  1 ,...  ,N. 

5.  Obtain  a  new  estimate  9rn+ 1  from  the  bootstrap  sample  {y™}  using  OLS.  Add  9m+1  into  the 
vector  0,  where  ©  is  a  vector  of  length  Mp  (M  is  the  number  of  bootstrap  samples)  which 
stores  the  bootstrap  estimates. 

6.  Set  m  =  m  +  1  and  repeat  steps  3-5. 

7.  Carry  out  the  above  iterative  process  M  times  where  M  is  large  (e.g.,  Af=1000),  resulting  in 
a  vector  0  of  length  Mp. 

8.  We  then  calculate  the  mean  and  standard  error  from  the  vector  0  using  the  formulae 

n  _  J_  sr^M  r\rn 

uboot  —  m  Z_^m=  1  u  ’ 

Cov (Oboot)  =  Em=  1  (0m  -  0boot)T(0m  -  Oboot),  (18) 

SEk(0boot)  =  \J  Cov  (Oboot)  kk- 

We  will  compare  the  optimal  design  methods  using  the  standard  errors  resulting  from  the  op¬ 
timal  time  points  each  method  proposes.  Since  there  are  different  ways  to  compute  the  standard 
errors  we  will  present  results  for  several  of  these  computational  methods. 
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4  The  Logistic  Growth  Example 


We  first  compare  the  optimal  design  methods  for  the  logistic  example  using  the  Monte  Carlo 
method  for  asymptotic  estimates  and  standard  errors. 


4.1  Logistic  Model 

The  Verhulst-Pearl  logistic  population  model  describes  a  population  that  grows  at  an  intrinsic 
growth  rate  until  it  reaches  its  carrying  capacity.  It  is  given  by  the  differential  equation: 

x(t)  =  rx(t)  ^1  -  ,  x(0)  =  xq, 


where  K  is  the  carrying  capacity  of  the  population,  r  is  the  intrinsic  growth  rate,  and  xo  is  the  initial 
population  size.  The  analytical  solution  to  the  differential  equation  above  is  given  by, 


x{t)  =  f(t,9  0) 


I< 

1  +  (K/x o  -  l)e_rf  ’ 


where  6q  =  (K.  r.  xq)  is  the  true  parameter  vector. 

Our  statistical  model  is  given  by 

Y(t)  =  f(t,0o) +£(t), 

where  we  choose  £  ~  Af(0,  Uq)  to  generate  simulated  data  (with  for  use  in  the  Monte  Carlo  calcula¬ 
tions).  A  realization  of  the  observation  process  is  given  by 


y(t)  =  f(t,90) +  e(t),  t  €  [0, T], 


4.2  Logistic  Results 

For  the  logistic  model,  we  use  SolvOpt  to  solve  for  the  optimal  mesh  for  each  of  the  optimal 
design  methods  (D-optimal,  FI-optimal  and  5'A-optimal) ,  using  the  second  constraint  (C 2)  on  the 
time  points:  t±  >  0,  tjy  <  T  and  1 ,  such  that  the  optimal  mesh  contains  0  and  T .  For  this 

example,  we  took  T  =  25  and  N  =  10  or  N  =  15.  Figures  1  and  3  contain  plots  of  the  resulting 
optimal  distribution  of  time  points  for  the  different  optimal  design  methods,  along  with  the  uniform 
mesh,  plotted  on  the  logistic  curve,  for  N  =  10  and  N  =  15,  respectively. 

These  optimal  design  methods  are  compared  based  on  their  average  Monte  Carlo  asymptotic 
estimates  and  standard  errors.  The  simulated  data  was  generated  assuming  the  true  parameter  values 
0o  =  (K,r,x o)  =  (17.5,0.8,0.1),  and  variance  <Tq  =  0.16.  The  average  estimates  and  standard  errors 
are  based  on  M  =  1000  Monte  Carlo  trials.  Since  we  obtain  histograms  of  estimates  and  standard 
errors  from  this  Monte  Carlo  analysis,  we  can  also  gain  information  for  comparison  from  the  median 
of  these  histograms  or  sampling  distributions.  Monte  Carlo  asymptotic  estimates  and  standard  errors 
were  also  computed  on  the  uniform  mesh.  We  report  the  average  and  median  estimates  and  standard 
errors  in  Tables  1  and  2  (. N  =  10,  and  N  =  15).  Histograms  for  the  Monte  Carlo  standard  errors  are 
given  in  Figs.  2  and  4  for  N  =  10  and  N  =  15  respectively. 

4.3  Discussion  of  Logistic  Results 

The  average  asymptotic  estimates  from  the  uniform  distribution  and  each  of  the  optimal  design 
methods  are  very  close  to  the  true  values,  9q.  For  N  =  10  (Table  1),  S’E'-optimal  has  the  closest 
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Table  1:  Average  and  Median  estimates  and  standard  errors  using  SolvOpt,  N  =  10,  M  =  1000,  and  do  =  (17.5,  0.7,  0.1). 
Optimization  with  constraint  implementation  (C 2). 


Parameter 

Method 

Average  Estimate 

Median  Estimate 

Average  SE 

Median  SE 

K 

Unif 

17.4978 

17.4954 

1.789  x  10'1 

1.789  x  10”1 

SE-  opt 

17.4985 

17.4995 

2.000  x  10”1 

2.000  x  10"1 

D-opt 

17.5066 

17.5024 

2.039  x  10'1 

2.038  x  10"1 

E-  opt 

17.4959 

17.4957 

1.512  x  10'1 

1.512  x  10”1 

r 

Unif 

0.7042 

0.6996 

5.020  x  10”2 

4.983  x  10-2 

SE-  opt 

0.7019 

0.7000 

3.473  x  10~2 

3.444  x  10~2 

H-opt 

0.7020 

0.7029 

3.821  x  10~2 

3.816  x  10~2 

E-  opt 

0.7139 

0.7033 

9.696  x  10”2 

9.090  x  10”2 

x0 

Unif 

0.1037 

0.0999 

3.730  x  10~2 

3.696  x  10~2 

SE- opt 

0.1018 

0.1002 

2.448  x  10-2 

2.432  x  10-2 

D-opt 

0.1025 

0.0982 

2.947  x  10”2 

2.859  x  10~2 

E-  opt 

0.1103 

0.0977 

6.417  x  10~2 

6.174  x  10"2 

Distribution  of  time  points 


Figure  1:  The  distribution  of  optimal  time  points  and  uniform  sampling  time  points  plotted  on  the  logistic  curve. 
Optimal  times  points  obtained  using  SolvOpt,  with  N  =  10,  and  the  optimal  design  methods  S' .E-optimality,  73- 
optimality,  and  75-optimality.  Optimization  with  constraint  implementation  (C 2). 


average  and  median  estimates,  followed  by  77-optimal  (for  r  and  xo)  and  E-optimal  (for  K ).  For 
N  =  15  (Table  2),  the  closest  average  estimate  of  K  came  from  E-optimal,  for  r  the  closest  average 
estimate  is  from  .SE-optimal  and  for  xo  it  was  D-optimal.  Comparing  the  average  and  median 
estimates,  we  see  that  for  all  cases  the  averages  and  medians  are  very  close,  indicating  that  the 
parameter  distributions  are  symmetric.  However,  for  both  N  =  10  and  N  =  15  the  averages 
were  slightly  larger  than  the  medians  for  r  and  xq  for  all  methods,  implying  that  those  parameter 
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Standard  Errors  for  K 


Standard  Errors  for  r 


(a) 


(b) 


Standard  Errors  for  xQ 


Figure  2:  Using  SolvOpt,  with  N  =  10,  a  comparison  of  optimal  design  methods  using  S' -E-optimality,  E-optimality, 
E-optimality,  with  a  uniform  sampling  time  points  in  terms  of  SEk  (panel  (a)),  SEr  (panel  (b)),  and  SEX o  (panel 
(c)).  Optimization  with  constraint  implementation  (C2). 


distributions  are  slightly  skewed  to  the  right  (see  Tables  1  and  2). 

Comparing  the  standard  errors  (Tables  1  and  2  and  Figs.  2  and  4):  For  K,  we  find  that  id- 
optimal  has  the  smallest  average  standard  errors,  then  the  uniform  grid,  then  S' id-optimal  when 
N  =  10  or  .D-optimal  when  N  =  15.  For  r  and  x'o,  Sid-optimal  has  the  smallest  average  standard 
errors,  followed  by  -D-optimal,  then  the  uniform  grid.  The  average  and  median  standard  errors  are 
very  close.  However  the  distribution  of  standard  errors  for  r  and  id-optimal  seem  to  be  slightly 
right-skewed. 

In  conclusion,  all  of  the  optimal  design  methods  produce  parameter  estimates  that  are  close  to 
the  true  value.  In  addition,  the  standard  error  estimates  are  similar  comparing  the  different  optimal 
design  methods.  Based  on  the  standard  errors,  id-optimal  is  more  favorable  for  the  accuracy  of  K , 
and  S' id-optimal  is  more  favorable  for  the  accuracy  of  r  and  xq  (followed  closely  by  -D-optimal). 
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Table  2:  Average  and  Median  estimates  and  standard  errors  using  SolvOpt,  N  =  15,  M  =  1000,  and  9q  =  (17.5,  0.7,  0.1). 
Optimization  with  constraint  implementation  (C2). 


Parameter 

Method 

Average  Estimate 

Median  Estimate 

Average  SE 

Median  SE 

K 

Unit 

17.5004 

17.5009 

1.467  x  10'1 

1.466  x  lO’1 

SE- opt 

17.4941 

17.4899 

1.633  x  10”1 

1.633  x  10”1 

H-opt 

17.5017 

17.4974 

1.553  x  10”1 

1.552  x  10"1 

E-  opt 

17.5006 

17.5006 

1.265  x  lO^1 

1.265  x  10”1 

r 

Unif 

0.7018 

0.6983 

4.118  x  10~2 

4.086  x  10-2 

SE- opt 

0.7008 

0.6993 

2.739  x  10~2 

2.721  x  10~2 

H-opt 

0.7022 

0.7016 

3.353  x  10~2 

3.353  x  10”2 

E-opt 

0.7056 

0.7004 

8.078  x  10”2 

7.799  x  10~2 

x0 

Unif 

0.1027 

0.1004 

3.040  x  10”2 

3.020  x  10”2 

SE- opt 

0.1016 

0.0997 

1.999  x  10~2 

1.977  x  10"2 

D-opt 

0.1014 

0.0992 

2.476  x  10”2 

2.440  x  10“2 

E-  opt 

0.1078 

0.0989 

4.920  x  10”2 

4.788  x  10”2 

Distribution  of  time  points 


Figure  3:  The  distribution  of  optimal  time  points  and  uniform  sampling  time  points  plotted  on  the  logistic  curve. 
Optimal  times  points  obtained  using  SolvOpt,  with  N  =  15,  and  the  optimal  design  methods  S' .E-optimality,  D- 
optimality,  and  E-optimality.  Optimization  with  constraint  implementation  (C2). 
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Standard  Errors  for  K 


Standard  Errors  for  r 


0  0.05  0.1  0.15  0.2  0.25 


(a) 


(b) 


Standard  Errors  for  xQ 


Figure  4:  Using  SolvOpt,  with  N  =  15,  a  comparison  of  optimal  design  methods  using  S' .E-optimality,  E-optimality, 
E-optimality,  with  a  uniform  sampling  time  points  in  terms  of  SEk  (panel  (a)),  SEr  (panel  (b)),  and  SEX o  (panel 
(c)).  Optimization  with  constraint  implementation  (C 2). 
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5  The  Harmonic  Oscillator  Model 


In  our  next  example,  we  consider  the  harmonic  oscillator,  also  known  as  the  spring-mass-dashpot 
model.  The  model  for  the  harmonic  oscillator  can  be  derived  using  Hooke’s  Law  and  mass-balance 
(see  [7])  and  is  given  by 


mx  +  cx  +  kx  =  0,  x(0)  =  x\,  x(0)  =  X2- 

Here,  m  is  mass,  c  is  damping,  and  k  is  the  spring  constant.  Dividing  through  by  m,  and  defining 
C  =  c/m  and  K  =  k/m ,  we  can  reduce  the  number  of  parameters. 

x  +  Cx  +  Kx  =  0,  x(0)  =  x'i ,  x(0)  =  X2- 

The  analytical  solution  for  the  position  at  time  t  can  be  obtained  and  is  given  by 

x(t)  =  e~at  (Ci  cos  bt  +  C2  sin  bt) , 

where  Ci  =  X2,  C2  =  (xi  +  0x2) /&,  a  =  |C,  and  b  =  y / K  —  \C2.  Substituting  in  Ci  and  C2,  we 
obtain, 

x(t)  =  x(t,  9q)  =  f(t,  80)  =  e~at  ^X2  cos  bt  +  1  — —  sin  bt^j  ,  for  0  <  t  <  T, 

where  for  our  considerations  the  true  parameter  vector  is  given  by  9q  =  (C,  K,  xi ,  X2)  =  (0.1, 0.2,  —1,  0.5) 
in  our  examples  here. 


5.1  Results  for  the  Oscillator  Model 

The  first  way  we  will  compare  these  optimal  design  methods,  given  that  we  know  8q  =  (C,  K,  xi ,  X2) 
=  (0.1, 0.2, —1,  0.5)  and  Uq  =  0.16,  is  to  simply  use  their  corresponding  standard  errors  from  the 
asymptotic  theory,  i.e. ,  the  values  of  SE(9q)  given  in  (14).  Recall  that  uncertainty  is  quantified 
by  constructing  confidence  intervals  using  parameter  estimates  with  the  asymptotic  standard  error. 
Since  our  main  focus  here  is  the  width  of  the  confidence  intervals,  we  can  forgo  the  obtaining  of  the 
parameter  estimates  themselves  which,  for  now,  we  tacitly  assume  may  be  similar  for  the  three  data 
sampling  distributions  we  investigate  here. 

The  optimal  time  points  for  each  of  the  three  optimal  design  methods  are  plotted  with  the  model 
for  different  T  and  N  under  the  first  constraint  implementation  (Cl)  in  Fig.  5,  the  second  constraint 
implementation  (C2)  in  Fig.  6,  the  third  constraint  implementation  (C 3)  in  Fig.  7,  and  the  last 
constraint  implementation  (C4)  in  Fig.  8.  The  standard  errors  (14)  from  the  asymptotic  theory 
corresponding  to  these  optimal  meshes  are  given  in  Table  3-6,  respectively  for  the  four  different 
constraints. 
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Optimal  mesh  with  N=5,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm.  Optimal  mesh  with  N=10,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm 


(c) 


(d) 


Figure  5:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  8q  =  (C,  K ,  xi,  *2), 
with  T  —  14.14  (one  period)  for  N  =  5  (panel  (a))  and  N  =  10  (panel  (c))  and  T  =  28.28  (two  periods)  for  IV  =  10 
(panel  (b))  and  N  =  20(panel  (d)).  Optimization  with  constraint  implementation  (Cl). 
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Table  3:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  do  =  (C,K,x  1,2:2),  optimization  with  constraint  implementation  (Cl). 


T 

N 

Method 

SE(C) 

SE(K) 

SE(x  1) 

SE{x2 ) 

14.14  (1-period) 

5 

S'E'-optimal 

7.603  x  10-2 

4.320  x  10-2 

2.869  x  10-1 

3.714  x  10'1 

5 

H-optimal 

8.244  x  10-2 

2.539  x  10-2 

2.551  x  10_1 

3.940  x  10”1 

5 

.E-optimal 

1.243  x  10"1 

2.508  x  10"2 

3.685  x  10"1 

3.815  x  10”1 

14.14  (1-period) 

10 

S'-E-optimal 

5.527  x  10-2 

2.519  x  10-2 

2.113  x  10-1 

2.716  x  lO^1 

10 

H-optimal 

5.963  x  10"2 

1.845  x  10"2 

1.949  x  10"1 

2.821  x  10"1 

10 

.E-optimal 

1.136  x  10_1 

4.187  x  10-2 

2.193  x  10_1 

2.272  x  10"1 

28.28  (2-periods) 

10 

S'-E-optimal 

4.049  x  10-2 

1.980  x  10-2 

2.604  x  10_1 

2.305  x  10'1 

10 

H-optimal 

3.919  x  10"2 

1.372  x  10"2 

1.936  x  10"1 

2.816  x  10"1 

10 

.E-optimal 

7.080  x  10-2 

2.343  x  10-2 

2.242  x  10-1 

2.274  x  10”1 

28.28  (2-periods) 

20 

S'-E-optimal 

2.438  x  10"2 

1.457  x  10"2 

1.517  x  10"1 

1.633  x  10”1 

20 

H-optimal 

3.177  x  10-2 

1.102  x  10-2 

1.609  x  10_1 

2.632  x  10”1 

20 

.E-optimal 

4.422  x  10"2 

1.608  x  10"2 

1.355  x  10"1 

1.385  x  10”1 

Table  4:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  do  =  (C,  K,  xi,  *2),  optimization  with  constraint  implementation  (C 2). 


T 

N 

Method 

SE(C) 

SE(I<) 

SE(x  1) 

SE{x2 ) 

14.14  (1-period) 

5 

S'-E-optimal 

7.900  x  10"2 

2.657  x  10"2 

2.852  x  10"1 

3.657  x  10”1 

5 

H-optimal 

8.251  x  10-2 

2.541  x  10-2 

2.561  x  10_1 

3.921  x  10”1 

5 

.E-optimal 

1.371  x  10_1 

2.900  x  10-2 

3.583  x  10"1 

3.736  x  10"1 

14.14  (1-period) 

10 

S'E'-optimal 

5.667  x  10"2 

2.484  x  10"2 

1.964  x  10"1 

2.310  x  10”1 

10 

H-optimal 

6.055  x  10-2 

1.648  x  10-2 

1.986  x  10_1 

2.822  x  10”1 

10 

.E-optimal 

8.507  x  10"2 

2.657  x  10"2 

2.211  x  10"1 

2.283  x  10"1 

28.28  (2-periods) 

10 

S'E'-optimal 

3.430  x  10-2 

2.149  x  10-2 

1.970  x  10_1 

2.274  x  10”1 

10 

H-optimal 

7.445  x  10-2 

1.711  x  10-2 

4.314  x  10-1 

3.919  x  10”1 

10 

.E-optimal 

8.826  x  10"2 

2.532  x  10"2 

2.132  x  10"1 

2.169  x  10"1 

28.28  (2-periods) 

20 

S'-E-optimal 

2.457  x  10-2 

1.500  x  10-2 

1.516  x  10_1 

1.784  x  10”1 

20 

H-optimal 

3.254  x  10"2 

1.166  x  10"2 

1.722  x  10"1 

2.867  x  10”1 

20 

.E-optimal 

5.135  x  10-2 

1.628  x  10-2 

1.451  x  10_1 

1.492  x  10"1 
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Optimal  mesh  with  N=5,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm.  Optimal  mesh  with  N=10,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm 


(c) 


(d) 


Figure  6:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  8q  =  (C,  K ,  xi,  *2), 
with  T  —  14.14  (one  period)  for  N  =  5  (panel  (a))  and  N  =  10  (panel  (c))  and  T  =  28.28  (two  periods)  for  IV  =  10 
(panel  (b))  and  N  =  20(panel  (d)).  Optimization  with  constraint  implementation  (C 2). 
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Optimal  mesh  with  N=5,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm.  Optimal  mesh  with  N=10,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm 


(c) 


(d) 


Figure  7:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  8q  =  (C,  K ,  xi,  *2), 
with  T  =  14.14  (one  period)  for  N  =  5  (panel  (a))  and  N  =  10  (panel  (c))  and  T  =  28.28  (two  periods)  for  IV  =  10 
(panel  (b))  and  N  =  20(panel  (d)).  Optimization  with  constraint  implementation  (C3). 
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Table  5:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  do  =  (C,K,x  1,2:2),  optimization  with  constraint  implementation  (C3). 


T 

N 

Method 

SE(C) 

SE(K) 

SE(x  1) 

SE{x2 ) 

14.14  (1-period) 

5 

S'E'-optimal 

7.900  x  10-2 

2.657  x  10-2 

2.852  x  10-1 

3.657  x  10'1 

5 

H-optimal 

9.483  x  10-2 

2.106  x  10-2 

2.675  x  10_1 

3.898  x  10”1 

5 

.E-optimal 

1.371  x  10"1 

2.900  x  10"2 

3.583  x  10"1 

3.736  x  10”1 

14.14  (1-period) 

10 

S'E'-optimal 

5.666  x  10-2 

2.484  x  10-2 

1.963  x  10_1 

2.309  x  10-1 

10 

D-optimal 

6.071  x  10-2 

1.656  x  10"2 

1.978  x  10"1 

2.828  x  10"1 

10 

.E-optimal 

1.125  x  10_1 

2.838  x  10-2 

2.532  x  10-1 

2.639  x  10"1 

28.28  (2-periods) 

10 

S'E'-optimal 

3.673  x  10-2 

2.399  x  10-2 

1.925  x  10_1 

2.000  x  10”1 

10 

D-optimal 

3.764  x  10-2 

1.373  x  10"2 

1.881  x  10"1 

2.812  x  10"1 

10 

.E-optimal 

7.949  x  10-2 

2.509  x  10-2 

2.154  x  10_1 

2.183  x  10”1 

28.28  (2-periods) 

20 

S'E'-optimal 

2.671  x  10"2 

1.812  x  10"2 

1.368  x  10"1 

1.413  x  10”1 

20 

D-optimal 

2.882  x  10-2 

1.057  x  10-2 

1.176  x  10-1 

1.959  x  10”1 

20 

.E-optimal 

6.467  x  10"2 

2.604  x  10"2 

1.361  x  10"1 

1.376  x  10"1 

Table  6:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  do  =  (C,  K,  xi,  *2),  optimization  with  constraint  implementation  (C4). 


T 

N 

Method 

SE{C) 

SE(I<) 

SE(x  1) 

SE{x2 ) 

14.14  (1-period) 

5 

S'E'-optimal 

7.900  x  10"2 

2.657  x  10"2 

2.852  x  10"1 

3.657  x  10”1 

5 

D-optimal 

8.249  x  10-2 

2.538  x  10-2 

2.553  x  10_1 

3.935  x  10”1 

5 

.E-optimal 

1.371  x  10-1 

2.900  x  10-2 

3.583  x  10"1 

3.736  x  10"1 

14.14  (1-period) 

10 

S'E'-optimal 

5.666  x  10"2 

2.484  x  10"2 

1.963  x  10"1 

2.309  x  10'1 

10 

D-optimal 

6.073  x  10-2 

1.657  x  10-2 

1.978  x  10-1 

2.828  x  10'1 

10 

.E-optimal 

1.125  x  10"1 

2.838  x  10"2 

2.532  x  10"1 

2.639  x  10"1 

28.28  (2-periods) 

10 

S'E'-optimal 

3.554  x  10-2 

2.395  x  10-2 

1.906  x  10_1 

2.000  x  10”1 

10 

D-optimal 

3.765  x  10-2 

1.373  x  10-2 

1.881  x  10-1 

2.812  x  10”1 

10 

.E-optimal 

7.948  x  10"2 

2.509  x  10"2 

2.154  x  10"1 

2.183  x  10"1 

28.28  (2-periods) 

20 

S'E'-optimal 

2.512  x  10-2 

1.698  x  10-2 

1.348  x  10_1 

1.510  x  10”1 

20 

D-optimal 

2.920  x  10"2 

1.035  x  10"2 

1.221  x  10"1 

1.788  x  10”1 

20 

.E-optimal 

6.095  x  10-2 

2.597  x  10-2 

1.300  x  10_1 

1.315  x  10"1 

20 


Optimal  mesh  with  N=5,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm.  Optimal  mesh  with  N=10,  and  0  =  (C,K,y  ,y  )  using  SolveOpt  algorithm 


(c) 


(d) 


Figure  8:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  8q  =  (C,  K ,  xi,  *2), 
with  T  —  14.14  (one  period)  for  N  =  5  (panel  (a))  and  N  =  10  (panel  (c))  and  T  =  28.28  (two  periods)  for  IV  =  10 
(panel  (b))  and  N  =  20(panel  (d)).  Optimization  with  constraint  implementation  (C4). 
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5.2  Discussion  for  the  Oscillator  Model 

The  constrained  optimization  algorithm,  SolvOpt,  was  chosen  over  MATLAB’s  fmincon  for 
comparisons  using  the  harmonic  oscillator  example  because  it  overall  resulted  in  more  well-behaved 
standard  errors  (real  and  finite  values),  and  fmincon  often  did  not. 

In  most  cases,  optimal  meshes  with  a  larger  number  of  points  were  nested  in  the  optimal  meshes 
with  a  reduced  the  number  of  points.  In  some  cases  for  T  =  28.28  (Figs  5  and  6)  doubling  the 
number  of  points  resulted  in  extra  points  being  dispersed  to  otherwise  empty  regions,  while  other 
points  were  nested  in  the  optimal  mesh  with  fewer  points.  Often  the  larger  number  of  points  in  the 
optimal  mesh  resulted  in  smaller  standard  errors. 

Examining  the  asymptotic  standard  errors  (Tables  3-6),  different  optimal  sampling  distributions 
produced  the  smallest  standard  errors  for  different  parameters,  with  no  optimal  design  method  having 
consistently  smaller  standard  errors.  For  C .  most  of  the  time  5'E-optimal  had  the  smallest  standard 
error,  then  ZD-optimal.  For  K ,  ZD-optimal  most  often  had  the  smallest  standard  error,  followed  by 
SE-optimal.  For  x\ ,  ID-optimal  had  the  smallest  standard  errors  in  most  cases.  For  X2,  either 
SE-optimal  or  E-optimal  had  the  smallest  standard  errors. 

The  standard  errors  from  the  different  optimal  design  methods  were  usually  on  the  same  order 
of  magnitude.  No  method  was  always  the  best  while  comparing  asymptotic  standard  errors,  though 
for  specific  parameters  some  optimal  sampling  distributions  were  favorable. 

Since  the  asymptotic  standard  errors  appear  explicitly  in  the  cost  function  we  are  minimizing  for 
S'E’-optimal  design,  it  may  not  be  fair  to  compare  these  methods  based  on  their  asymptotic  standard 
errors.  To  account  for  any  possible  bias  in  our  comparison,  we  will  compare  these  optimal  design 
methods  in  the  next  section  using  simulated  data  and  the  inverse  problem  to  estimate  parameters 
using  asymptotic  theory  and  bootstrapping.  In  these  computations,  we  will  compare  the  optimal 
design  methods  based  on  how  close  their  parameter  estimates  are  to  the  true  parameters,  and  the 
values  of  their  estimated  standard  errors  and  covariances. 

5.3  Results  for  the  Oscillator  Model  -  with  the  Inverse  Problem 

We  solve  the  inverse  problem  with  the  OLS  formulation  to  obtain  parameter  estimates  and 
standard  errors  from  both  asymptotic  theory  (15)  and  the  bootstrapping  method  (18).  We  create 
simulated  noisy  data  (in  agreement  with  our  statistical  model  (2))  corresponding  to  the  optimal  time 
meshes  using  true  values  8q  =  (C,  K,  x±,  X2)  =  (0.1,  0.2, —1, 0.5)  and  iid,  noise  with  £j  ~  A7(0,Oq). 
In  this  section  we  only  estimate  a  subset  of  the  parameters  9  =  (C,  K).  In  addition  to  the  estimates 
and  standard  errors,  we  also  report  the  estimated  Co v(C,K)  according  to  asymptotic  theory  (12) 
and  bootstrapping  (18).  For  comparison  purposes  we  also  present  these  results  for  a  uniform  grid 
using  the  same  T  and  N. 

The  optimal  time  points  for  each  of  the  three  optimal  design  methods  are  plotted  with  the  model 
for  T  =  14.14  and  T  =  28.28  for  N  =  15  under  the  first  constraint  implementation  (Cl)  in  Fig.  9, 
the  second  constraint  implementation  (C 2)  in  Fig.  10,  the  third  constraint  implementation  (C3)  in 
Fig.  11,  and  the  last  constraint  implementation  (C 4)  in  Fig.  12.  The  estimates,  standard  errors, 
and  covariance  between  parameters  as  estimated  from  the  asymptotic  theory  (15)  corresponding 
to  these  optimal  meshes  are  given  in  Table  7,  9,  11,  and  13,  respectively  for  the  four  different 
constraint  implementations.  The  estimates,  standard  errors,  and  covariance  between  parameters 
when  estimated  from  the  bootstrapping  method  (18)  corresponding  to  these  optimal  meshes  are 
given  in  Table  8,  10,  12,  and  14,  respectively  for  the  four  different  constraints.  In  each  of  the  tables 
are  also  results  on  the  uniform  grid  of  time  points  for  the  same  T  and  N.  Since  this  is  unaffected  by 
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constraints,  the  results  for  the  uniform  grid  are  repeated  in  the  tables. 


Optimal  mesh  with  N=15,  and  0  =  (C,K)  using  SolveOpt  algorithm. 


t 


Optimal  mesh  with  N=15,  and  0  =  (C,K)  using  SolveOpt  algorithm. 


t 


(a) 


(b) 


Figure  9:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  9q  =  ( C,K ), 
N  =  15,  with  T  =  14.14  (one  period)  (panel  (a))  and  T  =  28.28  (two  periods)  (panel  (b)).  Optimization  with  constraint 
implementation  (Cl). 


Table  7:  Estimates  and  standard  errors  from  the  asymptotic  theory  (15)  resulting  from  different  optimal  design  methods 
(as  well  as  for  the  uniform  mesh)  for  do  =  (C,  K)  =  (0.1,  0.2)  and  N  =  15,  optimization  with  constraint  implementation 
(Cl). 


T 

Method 

Casy 

SE(Casy) 

k 

lvasy 

SE(Kasy) 

Cov(CaSy,  Kasy) 

14.14 

SlE-optimal 

0.0865 

1.369  x  10"2 

0.1979 

1.165  x  10"2 

-3.597  x  10"5 

14.14 

P-optimal 

0.1112 

2.104  x  10-2 

0.2038 

8.974  x  10-3 

-1.027  x  10-4 

14.14 

C-optimal 

0.0592 

3.009  x  10"2 

0.1736 

1.285  x  10"2 

-9.801  x  10"5 

14.14 

Uniform 

0.1300 

3.529  x  10-2 

0.1938 

1.278  x  10-2 

-2.803  x  10"4 

28.28 

ST-optimal 

0.1111 

3.221  x  10"2 

0.2040 

2.827  x  10"2 

-3.391  x  10-4 

28.28 

P-optimal 

0.0705 

1.710  x  10"2 

0.1974 

7.444  x  10"3 

-6.045  x  10"5 

28.28 

C-optimal 

0.0843 

1.664  x  10-2 

0.1953 

1.381  x  10-2 

4.378  x  10”5 

28.28 

Uniform 

0.0854 

1.792  x  10"2 

0.2122 

7.326  x  10-3 

-6.219  x  10"5 

5.4  Discussion  of  Oscillator  Results  with  the  Inverse  Problem 

The  simulated  data  was  created  using  the  “true”  parameter  values  =  ( C,K )  =  (0.1,  0.2).  So 
we  can  compare  the  optimal  design  methods  based  on  how  close  the  parameter  estimates  are  as  well 
as  how  large  the  estimates  of  the  standard  errors  and  covariances  are. 

For  asymptotic  estimates: 

Comparing  optimal  design  methods  based  on  which  has  parameter  estimates  closest  to  the 
true  values,  there  is  no  method  that  is  always  the  best.  For  constraint  implementation  (Cl),  with 
T  =  14.14  (Table  7)  the  closest  parameter  estimates  result  from  either  S'U-optimal  or  O-optimal.  For 
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Table  8:  Estimates  and  standard  errors  from  the  bootstrap  method  (18)  resulting  from  different  optimal  design  methods 
(as  well  as  for  the  uniform  mesh)  for  do  =  ( C,K )  =  (0.1,  0.2),  M  =  1000  bootstraps  and  N  =  15,  optimization  with 
constraint  implementation  (Cl). 


T 

Method 

Cboot 

SE(Cboot) 

E-boot 

SE(kboot) 

Cov (Cbooti  E-boot) 

14.14 

S'C-optimal 

0.0871 

1.460  x  10-2 

0.1988 

1.092  x  10-2 

8.013  x  10"5 

14.14 

H-optimal 

0.1035 

1.565  x  10"2 

0.2025 

8.225  x  10"3 

-1.731  x  10"5 

14.14 

C-optimal 

0.0603 

2.861  x  10"2 

0.1743 

1.297  x  10"2 

6.383  x  10~5 

14.14 

Uniform 

0.1170 

2.469  x  10-2 

0.1989 

1.009  x  10-2 

-4.978  x  10"5 

28.28 

S'C-optimal 

0.0827 

2.486  x  10"2 

0.1991 

1.624  x  10"2 

1.567  x  10“4 

28.28 

H-optimal 

0.0705 

1.428  x  10-2 

0.1972 

7.177  x  10-3 

-5.355  x  10-6 

28.28 

U-optimal 

0.0843 

2.138  x  10"2 

0.2035 

2.598  x  10"2 

3.099  x  10"6 

28.28 

Uniform 

0.0837 

1.475  x  10"2 

0.2122 

6.350  x  10"3 

-4.436  x  10"6 

Optimal  mesh  with  N=15,  and  9  =  (C,K)  using  SolveOpt  algorithm. 


Optimal  mesh  with  N=15,  and  9  =  (C,K)  using  SolveOpt  algorithm. 


(a) 


(b) 


Figure  10:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  do  =  ( C,K ), 
N  =  15,  with  T  =  14.14  (one  period)  (panel  (a))  and  T  =  28.28  (two  periods)  (panel  (b)).  Optimization  with  constraint 
implementation  (C2). 


constraint  implementation  (C 2)  (Table  9),  either  D-optimal  or  U-optimal  had  the  closest  parameter 
estimates  to  the  true  values.  For  constraint  implementation  (C 3)  (Table  11),  either  H-optimal  or 
U-optimal  had  the  closest  parameter  estimate  for  C .  and  either  S'U-optimal  or  H-optimal  has  the 
closest  estimate  for  K .  For  constraint  implementation  (C4)  (Table  13),  SU-optimal  and  H-optimal 
had  the  closest  estimates  for  T  =  14.14,  and  U-optimal  had  the  closest  estimates  for  T  =  28.28. 

Comparing  the  optimal  design  methods  based  on  the  estimated  standard  errors  and  covariance 
between  parameters,  we  find  that  no  method  is  always  the  best.  For  constraint  implementation 
(Cl)  (Table  7),  when  T  =  14.14  5'C-optimal  had  the  smallest  standard  errors  and  covariance,  when 
T  =  28.28  either  C-optimal  or  D-optimal  had  the  smallest  standard  errors  and  covariances.  For 
constraint  implementation  (C2)  (Table  9),  the  smallest  standard  errors  and  covariances  came  from 
C-optimal  when  T  =  14.14  and  S'C-optimal  when  T  =  28.28,  and  followed  by  H-optimal  in  both 
cases.  For  constraint  implementation  (C3)  (Table  11),  the  smallest  standard  errors  and  covariances 
came  from  H-optimal  or  C-optimal  when  T  =  14.14,  and  from  H-optimal  followed  by  S'C-optimal 
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Table  9:  Estimates  and  standard  errors  from  the  asymptotic  theory  (15)  resulting  from  different  optimal  design  methods 
(as  well  as  for  the  uniform  mesh)  for  80  =  (C,  K)  =  (0.1,  0.2)  and  N  =  15,  optimization  with  constraint  implementation 
(C2). 


T 

Method 

Casy 

SE(Casy) 

K 

lvasy 

SE(Kasy) 

Cov(CaSy,  Kasy) 

14.14 

S'U-optimal 

0.0841 

2.852  x  10-2 

0.2314 

1.996  x  10-2 

-3.540  x  10-4 

14.14 

H-optimal 

0.0934 

2.635  x  10"2 

0.2054 

9.968  x  10"3 

-1.414  x  10~4 

14.14 

U-optimal 

0.1076 

2.220  x  10"2 

0.1952 

1.060  x  10"2 

-9.008  x  10~5 

14.14 

Uniform 

0.1300 

3.529  x  10-2 

0.1938 

1.278  x  10-2 

-2.803  x  10-4 

28.28 

S'U-optimal 

0.0649 

1.440  x  10"2 

0.1842 

7.006  x  10"3 

2.883  x  10~6 

28.28 

D-optimal 

0.1088 

1.888  x  10-2 

0.2086 

8.425  x  10-3 

-6.880  x  10-5 

28.28 

U-optimal 

0.1115 

2.397  x  10"2 

0.2046 

2.073  x  10"2 

-1.256  x  10~4 

28.28 

Uniform 

0.0854 

1.792  x  10"2 

0.2122 

7.326  x  10-3 

-6.219  x  10"5 

Table  10:  Estimates  and  standard  errors  from  the  bootstrap  method  (18)  resulting  from  different  optimal  design 
methods  (as  well  as  for  the  uniform  mesh)  for  6q  =  ( C,K )  =  (0.1,  0.2),  M  =  1000  bootstraps  and  N  =  15,  optimization 
with  constraint  implementation  (C 2). 


T 

Method 

C'boot. 

SE(Cboot) 

Eboot 

SE(Kboot) 

Co v(Cbooti  Eboot) 

14.14 

S'U-optimal 

0.0783 

2.076  x  10-2 

0.2320 

1.628  x  10-2 

4.751  x  10~5 

14.14 

H-optimal 

0.0976 

2.243  x  10"2 

0.2070 

9.921  x  10"3 

-4.040  x  10"5 

14.14 

U-optimal 

0.1031 

1.930  x  10-2 

0.1956 

9.636  x  10-3 

3.043  x  10~5 

14.14 

Uniform 

0.1170 

2.469  x  10"2 

0.1989 

1.009  x  10"2 

-4.978  x  10"5 

28.28 

S'U-optimal 

0.0576 

1.479  x  10-2 

0.1842 

6.057  x  10-3 

3.937  x  10”5 

28.28 

H-optimal 

0.1194 

1.694  x  10-2 

0.2105 

8.317  x  10-3 

4.750  x  10~6 

28.28 

U-optimal 

0.0947 

2.161  x  10"2 

0.2045 

1.927  x  10"2 

1.499  x  10~4 

28.28 

Uniform 

0.0837 

1.475  x  10-2 

0.2122 

6.350  x  10-3 

-4.436  x  10-6 

when  T  =  28.28.  For  constraint  implementation  (C 4)  when  T  =  14.14,  Fi-optimal  had  the  smallest 
standard  errors  and  covariances  followed  by  D-optimal. 

For  bootstrap  estimates: 

Comparing  optimal  design  methods  based  on  which  has  bootstrapping  parameter  estimates  clos¬ 
est  to  the  true  value,  again  no  method  is  always  the  best.  For  constraint  implementations  (Cl)  and 
(C 4)  (Tables  8  and  14),  when  T  =  14.14  either  S'Fh-optimal  or  D-optimal  have  the  closest  estimates. 
For  constraint  implementation  (C 2)  (Table  10),  either  D-optimal  or  F’-optimal  had  parameter  esti¬ 
mates  closest  to  the  true  values.  For  T  =  14.14  (Table  10),  the  parameter  estimate  for  K  was  in  fact 
closest  from  the  uniform  mesh,  followed  by  H-optimal.  For  constraint  implementation  (C 3)  (Table 
12),  when  T  =  14.14  either  D-optimal  or  U-optimal  had  the  closest  estimates.  For  cases  that  were 
not  reported,  there  was  no  method  that  was  consistently  better  in  terms  of  closeness  of  parameter 
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Optimal  mesh  with  N=15,  and  0  =  (C,K)  using  SolveOpt  algorithm.  Optimal  mesh  with  N=15,  and  0  =  (C,K)  using  SolveOpt  algorithm. 


(a)  (b) 


Figure  11:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  6 o  =  ( C,K ), 
N  =  15,  with  T  =  14.14  (one  period)  (panel  (a))  and  T  =  28.28  (two  periods)  (panel  (b)).  Optimization  with  constraint 
implementation  (C 3). 

Table  11:  Estimates  and  standard  errors  from  the  asymptotic  theory  (15)  resulting  from  different  optimal  design 
methods  (as  well  as  for  the  uniform  mesh)  for  do  =  ( C,K )  =  (0.1,  0.2)  and  N  =  15,  optimization  with  constraint 
implementation  (C 3). 


T 

Method 

Casy 

SE{Casy) 

K 

1Yasy 

SE(Kasy) 

Co ’v{Casy<)  -Kasy) 

14.14 

SC-optimal 

0.1238 

2.515  x  10"2 

0.2028 

2.302  x  10"2 

-1.640  x  10"4 

14.14 

C-optimal 

0.0970 

2.061  x  10-2 

0.1997 

7.973  x  10-3 

-9.382  x  10-5 

14.14 

C-optimal 

0.1156 

2.204  x  10"2 

0.1953 

2.055  x  10"2 

6.635  x  10”5 

14.14 

Uniform 

0.1300 

3.529  x  10-2 

0.1938 

1.278  x  10-2 

-2.803  x  10-4 

28.28 

SC-optimal 

0.0899 

1.617  x  10-2 

0.2015 

1.368  x  10-2 

-5.288  x  10-5 

28.28 

C-optimal 

0.0966 

1.540  x  10"2 

0.2084 

6.787  x  10"3 

-4.422  x  10~5 

28.28 

C-optimal 

0.1029 

1.705  x  10-2 

0.2098 

2.111  x  10-2 

-1.575  x  10-4 

28.28 

Uniform 

0.0854 

1.792  x  10"2 

0.2122 

7.326  x  10-3 

-6.219  x  10"5 

estimates  to  the  true  values. 

Comparing  optimal  design  methods  based  on  which  method  produces  the  smallest  bootstrap¬ 
ping  estimated  standard  errors  and  parameter  estimates,  no  method  is  consistently  favorable.  For 
constraint  implementation  (Cl)  (Table  8),  C-optimal  has  the  smallest  standard  errors  and  covari¬ 
ances.  For  constraint  implementation  (C 2)  (Table  10),  when  T  =  14.14  the  smallest  standard  errors 
and  covariances  come  from  C-optimal,  when  T  =  28.28  either  SC-optimal  or  the  uniform  grim  had 
the  smallest  standard  errors  and  covariances,  followed  by  C-optimal.  For  constraint  implementa¬ 
tion  (C3)  (Table  12),  the  smallest  standard  errors  and  covariances  are  from  C-optimal,  followed  by 
SC-optimal.  For  constraint  implementation  (C4)  (Table  14),  when  T  =  14.14  the  smallest  stan¬ 
dard  errors  and  covariances  come  from  C-optimal  followed  by  C-optimal,  when  T  =  28.28  either 
SC-optimal  or  C-optimal  were  the  smallest. 

In  conclusion,  all  of  the  optimal  design  methods  are  favorable  under  specific  conditions.  In 
many  of  the  cases  the  parameter  estimates,  standard  errors,  and  covariances  are  on  the  same  order 


26 


Table  12:  Estimates  and  standard  errors  from  the  bootstrap  method  (18)  resulting  from  different  optimal  design 
methods  (as  well  as  for  the  uniform  mesh)  for  do  =  ( C,K )  =  (0.1,  0.2),  M  =  1000  bootstraps  and  N  =  15,  optimization 
with  constraint  implementation  (C3). 


T 

Method 

14.14 

S'if-optimal 

0.1204 

2.652  x  10-2 

0.2047 

2.186  x  10-2 

3.199  x  10"4 

14.14 

D-optimal 

0.0919 

1.574  x  10"2 

0.1978 

7.301  x  10-3 

-2.363  x  10"5 

14.14 

U-optimal 

0.1069 

2.756  x  10"2 

0.1978 

1.967  x  10"2 

3.763  x  10~4 

14.14 

Uniform 

0.1170 

2.469  x  10-2 

0.1989 

1.009  x  10-2 

-4.978  x  10"5 

28.28 

S'U-optimal 

0.0870 

1.753  x  10"2 

0.2028 

1.542  x  10"2 

8.280  x  10~5 

28.28 

D-optimal 

0.0906 

1.133  x  10-2 

0.2030 

5.540  x  10-3 

5.491  x  10“6 

28.28 

U-optimal 

0.0926 

1.783  x  10"2 

0.2113 

2.202  x  10"2 

4.422  x  10~5 

28.28 

Uniform 

0.0837 

1.475  x  10"2 

0.2122 

6.350  x  10"3 

-4.436  x  10"6 

Optimal  mesh  with  N=15,  and  9  =  (C,K)  using  SolveOpt  algorithm.  Optimal  mesh  with  ISM  5,  and  9  =  (C,K)  using  SolveOpt  algorithm. 


(a)  (b) 

Figure  12:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  do  =  ( C,K ), 
N  =  15,  with  T  =  14.14  (one  period)  (panel  (a))  and  T  =  28.28  (two  periods)  (panel  (b)).  Optimization  with  constraint 
implementation  (C4). 


of  magnitude  resulting  from  different  optimal  design  criteria. 
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Table  13:  Estimates  and  standard  errors  from  the  asymptotic  theory  (15)  resulting  from  different  optimal  design 
methods  (as  well  as  for  the  uniform  mesh)  for  #o  =  ( C,K )  =  (0.1,  0.2)  and  N  =  15,  optimization  with  constraint 
implementation  (C4). 


T 

Method 

SE{CaSy) 

SE(Kasy) 

Cov(CaSy,  -Kasy) 

14.14 

S'U-optimal 

0.0837 

2.312  x  10"2 

0.2136 

2.087  x  10"2 

-2.385  x  10"4 

14.14 

D-optimal 

0.0771 

2.206  x  10-2 

0.1827 

8.461  x  10-3 

-1.074  x  10”4 

14.14 

£,-optimal 

0.0343 

1.387  x  10-2 

0.1719 

7.166  x  10-3 

2.825  x  10~5 

14.14 

Uniform 

0.1300 

3.529  x  10"2 

0.1938 

1.278  x  10"2 

-2.803  x  10"4 

28.28 

S'U-optimal 

0.0908 

1.473  x  10-2 

0.2206 

1.330  x  10-2 

-1.480  x  10“4 

28.28 

ZUoptimal 

0.1160 

2.358  x  10"2 

0.1875 

9.149  x  10"3 

-9.691  x  10"5 

28.28 

U-optimal 

0.0964 

1.218  x  10-2 

0.2070 

1.445  x  10-2 

-4.984  x  10~5 

28.28 

Uniform 

0.0854 

1.792  x  10-2 

0.2122 

7.326  x  10-3 

-6.219  x  10-5 

Table  14:  Estimates  and  standard  errors  from  the  bootstrap  method  (18)  resulting  from  different  optimal  design 
methods  (as  well  as  for  the  uniform  mesh)  for  do  =  ( C,K )  =  (0.1,  0.2),  M  =  1000  bootstraps  and  N  =  15,  optimization 
with  constraint  implementation  (C4). 


T 

Method 

14.14 

S'U-optimal 

0.0856 

2.185  x  10-2 

0.2184 

2.027  x  10-2 

1.904  x  10”4 

14.14 

D-optimal 

0.0658 

1.611  x  10"2 

0.1822 

7.297  x  10"3 

-2.611  x  10"5 

14.14 

U-optimal 

0.0334 

1.769  x  10-2 

0.1729 

6.841  x  10-3 

7.838  x  10-5 

14.14 

Uniform 

0.1170 

2.469  x  10-2 

0.1989 

1.009  x  10-2 

-4.978  x  10"5 

28.28 

S'U-optimal 

0.0835 

8.868  x  10"3 

0.2111 

7.826  x  10"3 

2.677  x  10~5 

28.28 

D-optimal 

0.1265 

1.872  x  10-2 

0.1986 

9.266  x  10-3 

-1.156  x  10"5 

28.28 

U-optimal 

0.0963 

1.594  x  10"2 

0.2195 

2.443  x  10"2 

-1.188  x  10~4 

28.28 

Uniform 

0.0837 

1.475  x  10-2 

0.2122 

6.350  x  10-3 

-4.436  x  10"6 
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6  A  Simple  Glucose  Regulation  Model 


Next  we  will  consider  a  well-known  model  for  the  intervenous  glucose  tolerance  test  (IVGTT). 
This  model  is  referred  to  as  the  minimal  model  in  the  literature  [10,  13,  21].  Prior  to  the  IVGTT 
the  patient  is  asked  to  fast.  When  the  patient  comes  in  for  the  IVGTT,  measurements  of  their 
baseline  glucose  and  insulin  concentrations,  Gb  and  Ib,  respectively,  are  first  taken.  The  IVGTT 
procedure  consists  of  injecting  a  bolus  resulting  in  an  initial  concentration  po  of  glucose  into  the 
blood,  and  measuring  the  glucose  and  insulin  concentrations  in  the  blood  at  various  time  points 
after  the  injection. 

The  body  carefully  regulates  the  glucose  concentration  in  the  blood  within  a  narrow  range. 
Extremely  high  blood  glucose  concentration  is  referred  to  as  hyperglycemia,  whereas  hypoglycemia 
results  when  the  blood  glucose  concentration  is  too  low.  The  IVGTT  initially  brings  the  blood 
glucose  concentration  to  hyperglycemic  levels.  In  normal  healthy  patients,  the  high  level  of  glucose 
in  the  blood  signals  the  beta  cells  of  the  pancreas  to  secrete  insulin.  Insulin  helps  the  fat  and  muscle 
cells  to  uptake  glucose  from  the  blood,  either  for  fuel  or  for  storage  as  glycogen.  When  the  blood 
glucose  concentration  is  too  low,  the  pancreas  secretes  glucagon  which  releases  glucose  stored  in  the 
liver  into  the  blood.  Glucagon  is  another  dynamic  variable  [3]  during  the  IVGTT.  Though  glucagon 
is  not  included  in  this  model,  it  is  acknowledged  that  the  liver  can  regulate  glucose  independently 
from  insulin  through  glucagon. 


6.1  Model 

The  minimal  model  is  given  by  the  following  system  of  ordinary  differential  equations  (see 
[10,  13,  21]  for  details): 

G(t)  =  -pi  ( G(t )  -  Gb)  -  X(t)G(t),  G( 0)  =  po,  (19) 

X(t)  =  -p2x(t)  +  p3  (/(*)  -  h) ,  x(0)  =  0,  (20) 

i(t)  =  p4tmax(0,G(i)  -  p5)  -p6(l(t)  -  /&),  I(0)=p7  +  Ib,  (21) 


where  G(t)  is  the  glucose  concentration  (in  mg/dl)  in  plasma  at  time  t,  I(t )  is  the  insulin  concentra¬ 
tion  (in  /iU/ml)  in  plasma  at  time  t  and  X(t)  represents  insulin-dependent  glucose  uptake  activity 
(proportional  to  a  remote  insulin  compartment)  in  units  1/min. 

We  use  the  following  approximate  max  function  in  equation  (21)  since  it  is  continuously  differ¬ 
entiable: 


maxfunci  (v) 


v  for  v  >  eo, 

<  0  for  v  <  — eo, 

T~  (w  +  eo)2  for  v  €  [— e0,e0], 
4eo 


where  eo  >  0  is  chosen  sufficiently  small  (for  example,  eo  =  ICE5). 

An  interpretation  of  the  parameters  is  given  in  Table  15. 

In  the  following  we  will  describe  the  model  and  its  underlying  assumptions. 


Equation  (19)  (Glucose  concentration  in  plasma) 

At  t  =  0  a  bolus  of  glucose  is  injected  such  that  the  initial  glucose  concentration  in  the  blood 
is  po-  The  first  term  represents  hepatic  glucose  balance,  which  occurs  independent  of  insulin  level. 
The  second  term  is  the  loss  of  glucose  due  to  insulin-dependent  uptake  by  peripheral  tissues. 
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Table  15:  Description  of  model  parameters  and  typical  values. 


e 

Description 

value 

Gb 

basal  pre-injection  level  of  glucose 

83.7  mg/dl 

h 

basal  pre-injection  level  of  insulin 

11  /iU/ml 

Po 

the  theoretical  glucose  concentration  in  plasma  at 

time  t  =  0 

279  mg/dl 

Pi 

the  rate  constant  of  insulin-independent  glucose 
uptake  in  muscles,  and  adipose  tissue 

2.6  x  10~2  min-1 

P2 

the  rate  constant  for  decrease  in  tissue  glucose 
uptake  ability 

0.025  min-1 

P3 

the  rate  constant  for  the  insulin-dependent  increase 
in  glucose  uptake  ability  in  tissue  per  unit  of  insulin 
concentration  above  Ib 

1.25  x  10-5  min_2(/iU/ml)“1 

Pi 

the  rate  constant  for  insulin  secretion  by  the 
pancreatic  /3-cells  after  the  glucose  injection  and 
with  glucose  concentration  above  pb 

4.1  x  lO-3  (/iU/ml)  min-2(mg/dl)-1 

Po 

the  threshold  value  of  glucose  in  plasma  above 
which  the  pancreatic  /3-cells  secrete  insulin 

83.7  mg/dl 

P6 

the  first  order  decay  rate  for  insulin  in  plasma 

0.27  min-1 

P7 

pi  +  Ib  is  the  theoretical  insulin  concentration  in 
plasma  at  time  t  =  0 

352.7  /xU/ml 

Equation  (20)  (Insulin-dependent  glucose  uptake  activity) 

At  t  =  0  there  is  no  glucose  uptake  activity.  Spontaneously  tissue  loses  the  ability  to  uptake 
glucose,  even  in  the  presence  of  insulin.  Glucose  uptake  activity  increases  proportionally  to  the 
amount  by  which  insulin  concentration  is  greater  than  baseline  insulin  concentration. 


Equation  (21)  (Insulin  concentration  in  the  plasma) 

At  t  =  0  the  initial  insulin  concentration  is  at  some  level  over  baseline,  given  by  pi  +  I b. 
The  increase  in  insulin  concentration  is  proportional  to  the  amount  by  which  glucose  concentration 
exceeds  some  threshold,  pb,  and  the  amount  of  time  that  has  elapsed  since  the  glucose  injection. 
There  is  a  loss  of  insulin  to  degradation  in  the  plasma.  The  pancreas  secretes  low  levels  of  insulin, 
even  in  hypoglycemic  conditions,  to  maintain  insulin  concentration  at  or  above  baseline  Ib. 

The  analysis  of  this  model  found  in  [10,  21]  gives  a  metabolic  portrait  for  the  first  phase  sensitiv¬ 
ity  to  glucose  ((pi)  (corresponding  to  initial  secretion  of  insulin),  the  second  phase  glucose  sensitivity 
( Sq )  (corresponding  to  a  secondary  phase  of  insulin  secretion),  and  the  insulin  sensitivity  index  (Si). 
The  metabolic  portrait  which  is  given  by 


Cl  tr  O  c<  JL 

ji  =  — ,  Sq  =  p  i,  (pi 

P2 


h 


P6  (po  -  Gb)  ’ 


(22) 
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where  Imax  is  the  maximal  value  of  insulin  concentration  in  plasma. 

Bergman  et  al.  [9]  suggest  the  use  of  this  model  in  the  clinical  IVGTT  setting.  Parameters 
from  the  model  are  estimated  using  patient-specific  data.  The  parameter  estimates  are  then  used 
in  the  metabolic  portrait  for  diabetes  diagnosis  purpose  for  that  patient.  This  process  was  made 
readily  available  to  clinicians  in  the  computer  software  MINMOD  [18].  Since  the  estimation  of  these 
parameters  plays  such  a  crucial  role  in  the  diagnosis,  it  appears  that  optimal  design  methods  would 
be  of  great  assistance.  Data  sampled  at  the  optimal  time  points  would  result  in  the  most  accurate 
metabolic  portrait  produced  by  this  mathematical  model. 

Next  we  will  describe  the  corresponding  statistical  model  for  this  system  involving  vector  ob¬ 
servations.  We  obtain  numerical  solutions  using  MATLAB’s  ode45  since  there  does  not  exist  an 
analytical  solution  to  this  system  of  differential  equations.  Let  z(t,  #o)  =  (G(t,  Oq ),X(t,  do),  I(t,  do))T 
represent  our  model  solution.  Since  we  can  observe  realizations  of  G(t,do)  and  I(t,9o),  but  not 
X(t,  9q),  our  observation  process  is  given  by 

m  =  f(t,e0)  =  (G(t,e0),i(t,e0))T. 

Our  statistical  model  is  given  by  the  stochastic  process 

Y(t)  =  f(t,6Q)  +  £(t), 

where  £{t)  is  a  noisy  vector  random  process.  We  assume  the  following  about  the  vector  random 
variable  £{t): 

E(£(t))  =  0,  ie[0,T], 

Var  £(t)  =  diag(o'o)G,  Oqj),  i€[0,T], 

Cov(£i(t)£i(s))  =  alG8(t-s),  t,se[0,T], 

Co v(£2(t)£2(s))  =  croj6(t-s),  t,s  £  [0, T], 

Cov(£i(t)£2(s))  =  0,  t,se[0,T], 

We  assume  constant  variance,  <Tq  g  =  25  and  <Tq  1  =  4.  A  realization  of  the  observation  process  is 
given  by 

y(t)  =  f{t,d0)  +  e(t),  t  <E  [0,  T], 
where  the  measurement  error  s(t)  is  a  realization  of  £{t). 

6.2  Methods 

Though  the  vector  methodology  is  similar  to  that  in  the  scaler  case,  for  completeness  we  outline 
it  here  for  a  system  of  differential  equations  such  as  the  simple  glucose  regulation  model. 

We  begin  by  finding  the  optimal  discrete  sampling  distribution  of  time  points  r  =  for 

a  fixed  number  of  points,  N,  and  a  fixed  final  time,  T,  using  either  S'A-optimal,  D-optimal,  or  Li- 
optimal.  These  three  optimal  design  methods  are  then  compared  based  on  the  asymptotic  standard 
errors  formulae  for  parameters  using  these  sampling  times. 

More  specifically,  once  we  have  an  optimal  distribution  of  time  points  we  will  obtain  data  or 
simulated  data,  {yi}f=n  a  realization  of  the  random  process  =  {(Gi,  Ii)T}^=1  given  by 

Yi  =  f(u,e  0)  +  £h 
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corresponding  to  the  optimal  time  points,  r  =  where  =  £{U). 

Define  V0  =  Var (Si)  =  diag(<ToiG,  Oq,/)- 

A  subset  of  the  parameters  is  estimated  by  inverse  problem  methodology  [4].  Since  the  variance 
is  assumed  to  be  constant,  the  inverse  problem  is  formulated  using  ordinary  least  squares  (OLS). 
The  OLS  estimator  for  a  vector  system  is  defined  by 

N 

©ols  =  ©ols  =  argmin^[y3  -  f(tj,  -  f(tj,6)\. 

3= 1 


For  a  given  realization  {i/j},  the  OLS  estimate  9qls  is  defined  as 


N 


^OLS  =  #ols  =  argmin^^  -  f(tj,9)]V0  1[yj  -  f(tj,6)]. 


3= 1 


The  definition  of  variance  gives 

Do  =  diag  E  (  1  J2[Yj  -  f(tj}9 0)}[Yj  -  f(tj}0 o)]3 


N 

j=1 

In  the  case  that  Vo  is  unknown  an  unbiased  estimate  can  be  obtained  from  the  realizations  {yi}f= i 
and  9  by 

(  1  N  _  .  . 

Eo  ~  V  =  diag  Y,Wi  ~  f(tj,  0)]Wj  ~  f(tjJ)]T 

\  p  j=i 

which  is  solved  simultaneously  (in  an  iterative  procedure  -  see  [4])  with  normal  equations  for  the 
estimate  9  =  9 olSj  where  p  is  the  number  of  parameters  being  estimated. 

To  compute  the  standard  errors  of  the  estimated  parameters,  we  first  must  compute  the  2  x  p 
sensitivity  matrices  Dj(9 )  =  (9)  which  are  given  by 


j 

dfi  {tj  ,9) 

dfi{tj,9)  \ 

Dj=( 

ddi 

d92 

99p  \ 

df2(tj,0) 

9  hit j, 9) 

dMtifi)  1  ’ 

\ 

\  d<h 

962 

99p  J 

for  j  =  1, . . . 

,  N.  For  this  system  we  can  rewrite  Dj  in 

terms  of  ( G(t.j,9)1I(t.j,9))T 

(since 

0),f2(tj,6))T  =  (G(tj,6),I(tj, 9)) 

\T).  We  have 

/  dG(th0) 

BGftjfi) 

9G(tj,8)  \ 

Dj  = 

|  0(h 

982 

99p  \ 

\  dl(tj,0) 

dl(tj,6) 

9i(tj,e) 

\  dOi 

dd2 

99p  J 

Y'AT 

^0 


The  true  covariance  matrix  is  approximately  (asymptotically  as  N  — >•  oo)  given  by 

-l 

\ 

YDJ(0o)vo~lDi(0o) 

0=1 

When  the  true  values,  9q  and  Vo,  are  unknown,  the  covariance  matrix  is  estimated  by 

-l 

\ 

Y  DJ (^OLsW^1  Dj {OoLs) 

0=1 

The  corresponding  FIM,  asymptotic  standard  errors  and  asymptotic  distribution  are  again  given 
by  (13),  (14),  (15),  and  (16),  respectively. 


«  tN  = 
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6.2.1  The  Bootstrap  Method  for  a  system 


The  bootstrap  method  for  a  system  of  differential  equations  is  the  same  as  described  in  the 
previous  section,  except  that  each  state  variable  has  its  own  residuals  that  must  be  separately  sam¬ 
pled  with  replacement.  The  first  four  steps  of  the  bootstrap  algorithm  of  Section  3.3  modified  for  a 
system  with  vector  observations  is  outlined  here  for  completeness. 


1.  First  estimate  9 0  from  the  entire  sample,  using  OLS. 

2.  Using  this  estimate  define  the  standardized  residuals: 


rG,j 


/2M0)) 


for  j  =  1 .... , N.  Then  {rcq,. . .  ,rG,N},{i"i, i,-  ■  ■  ivi,n}  are  realizations  of  rid  random  variables 
from  the  empirical  distribution  Jqy,  and  p  for  the  number  of  parameters. 


Set  m  =  0. 


3.  Create  a  two  different  bootstrap  samples  of  size  N  using  random  sampling  with  replacement 
from  the  data  (realizations)  {rcq,. . .  ,vg,n}  and  {f/q,. . .  ,77^}  to  form  the  bootstrap  samples 
irG,i’  •  •  •  >  vg,n}  and  {ry)l7 . . . ,  t™n}- 

4.  Create  bootstrap  sample  points 


j/iWiMV’-cj. 

yl,  =  h{t„t)")  +  rtP 

where  j  =  1,. . .  ,N. 

5.  Steps  5-8  are  the  same  as  those  of  the  algorithm  for  scalar  observations  in  Section  3.3. 

We  compute  the  optimal  time  mesh  using  S'E'-optimality,  D-optimality,  and  U-optimality,  as 
defined  in  the  previous  section,  for  a  subset  of  the  parameters  9  =  (pi,P2,P3,Pa),  and  a  fixed  number 
of  time  points  (N  =  30)  and  a  final  time  of  T  =  150  minutes.  We  remark  that  a  subset  of  parameters 
was  chosen  to  avoid  an  ill-conditioned  FIM.  The  subset  of  parameters  was  chosen  based  on  the 
traditional  sensitivity  functions.  The  glucose  and  insulin  model  solutions  were  most  sensitive  to 
0  =  ( P\iP2iPz-,Pa )•  The  approximate  asymptotic  standard  errors  (14)  for  6  =  (pi,P2,P3,P4)  were 
computed  on  the  optimal  mesh  corresponding  to  an  optimal  design  method. 

The  optimal  design  methods  were  implemented  using  the  constrained  minimization  algorithm 
SolvOpt.  The  variations  on  the  constraint  employed  were  the  same  as  in  the  previous  section. 
We  compare  SU-optimal,  H-optimal  and  U-optimal  design  methods  based  on  these  approximate 
asymptotic  standard  errors. 
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6.3  Results  for  the  Glucose  Regulation  Model 

The  optimal  time  points  (found  using  the  SolvOpt  algorithm)  for  each  of  the  three  optimal 
design  methods  are  plotted  with  the  model  for  T  =  150  minutes  and  N  =  30  under  the  first 
constraint  implementation  (Cl)  in  Fig.  13,  the  second  constraint  implementation  (C 2)  in  Fig.  14, 
the  third  constraint  implementation  (C3)  in  Fig.  15,  and  the  last  constraint  implementation  (C4) 
in  Fig.  16.  The  standard  errors  (14)  from  the  asymptotic  theory  corresponding  to  these  optimal 
meshes  are  given  in  Table  16-19,  respectively  for  the  four  different  constraint  implementations. 

Note  that  for  constraint  implementations  (C 2)  and  (C4)  initializing  SolvOpt  with  the  uniform 
mesh  resulted  in  a  terminal  error  for  C-optimal,  stating  that  the  gradient  at  the  starting  point  was 
zero.  In  these  cases  other  initial  mesh  points  were  chosen  such  that  C-optimal’s  initial  gradient  was 
non-zero,  and  optimization  could  be  achieved.  To  be  consistent,  all  three  methods  were  initialized  by 
the  same  non-uniform  mesh.  For  (C 2)  the  initial  mesh  was  r°  =  {0, . . . ,  0, 10,  37, 150, . . .  ,  150},  and 
for  (C4)  it  was  r°  =  {5, 15, 19,  21,  24,  26, 42,  59,  63,  73,  82,  95,  98,  98, 102,  111,  114, 119, 120, 122, 127, 
136, 137, 137, 140, 144, 144, 144, 145, 146}.  Optimal  design  methods  are  guaranteed  to  converge  in  a 
local  sense. 


Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm.  Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm. 


(a)  (b) 


Figure  13:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  9q  = 
(pi,P2,P3,P4),  with  T  =  150  for  N  =  30.  Optimal  time  points  with  the  Glucose  model  in  panel  (a)  and  with  the 
Insulin  model  in  panel  (b).  Optimization,  using  SolvOpt,  with  constraint  implementation  (Gl). 


Table  16:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  9q  =  (pi,p2,P3,P4),  optimization,  using  SolvOpt,  with  constraint  implementation  (Gl). 


Method 

SE(p  r) 

SE(p  2) 

SE(p3) 

SE(Pi) 

S'C-optimal 

4.173  x  10”3 

6.501  x  10”3 

3.100  x  10-6 

2.959  x  10-4 

C-optimal 

8.411  x  10"3 

1.236  x  10~2 

6.133  x  10-6 

1.714  x  10"4 

C-optimal 

4.381  x  10”3 

6.520  x  10”3 

3.182  x  10-6 

4.941  x  10-4 
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Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm.  Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm. 


t  t 


(a)  (b) 

Figure  14:  Plot  of  model  with  optimal  time  points  resulting  from  different,  optimal  design  methods  for  do  = 
(pi,P2,P3,P4),  with  T  =  150  for  N  =  30.  Optimal  time  points  with  the  Glucose  model  in  panel  (a)  and  with  the 
Insulin  model  in  panel  (b).  Optimization,  using  SolvOpt,  with  constraint  implementation  (C 2). 

Table  17:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  do  =  (pi,P2,P3,P4),  optimization,  using  SolvOpt,  with  constraint  implementation  (C 2). 


Method 

SE(p  1) 

SC(p2) 

SE(p3) 

sc(p4) 

SC-optimal 

4.019  x  10”3 

6.451  x  10”3 

3.088  x  10-6 

3.452  x  10-4 

C-optimal 

8.322  x  10~3 

1.103  x  10~2 

6.230  x  10"6 

2.748  x  10"4 

C-optimal 

3.882  x  10~3 

6.284  x  10~3 

3.063  x  10-6 

5.390  x  10"4 

Table  18:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  do  =  (pi,P2,P3,P4),  optimization,  using  SolvOpt,  with  constraint  implementation  (C3). 


Method 

SE(p  1) 

SE(p2 ) 

SE(p  3) 

sc(p4) 

SC-optimal 

4.205  x  10~3 

6.535  x  10”3 

3.151  x  10-6 

3.041  x  10-4 

C-optimal 

7.434  x  10”3 

1.517  x  10~2 

6.171  x  10"6 

1.181  x  10"4 

C-optimal 

7.528  x  10~3 

1.123  x  10~2 

5.509  x  10-6 

1.833  x  10"4 

6.4  Discussion  for  the  Glucose  Regulation  Model 

Comparing  the  optimal  design  methods  using  approximate  asymptotic  standard  errors,  we  find 
that  the  optimal  design  methods  that  are  best  for  (pi,P2,P3)  are  different  than  the  ones  best  for 
the  standard  error  of  p\.  For  constraint  implementation  (Cl)  (Table  16),  SC-optimal  followed  by 
C-optimal  had  the  smallest  standard  errors  for  (p\,p2,p:i),  and  C-optimal  followed  by  SC-optimal 
had  the  smallest  standard  errors  for  p 4.  For  constraint  implementation  (C 2)  (Table  17),  the  smallest 
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Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm.  Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm. 


(a)  (b) 


Figure  15:  Plot  of  model  with  optimal  time  points  resulting  from  different,  optimal  design  methods  for  do  = 
(pi,P2,P3,P4),  with  T  =  150  for  N  =  30.  Optimal  time  points  with  the  Glucose  model  in  panel  (a)  and  with  the 
Insulin  model  in  panel  (b).  Optimization,  using  SolvOpt,  with  constraint  implementation  (C3). 


Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm.  Optimal  mesh  with  N=30,  and  T=150  using  SolveOpt  algorithm. 


(a)  (b) 


Figure  16:  Plot  of  model  with  optimal  time  points  resulting  from  different  optimal  design  methods  for  do  = 
(pi,P2,P3,P4),  with  T  =  150  for  N  =  30.  Optimal  time  points  with  the  Glucose  model  in  panel  (a)  and  with  the 
Insulin  model  in  panel  (b).  Optimization,  using  SolvOpt,  with  constraint  implementation  (C4). 

standard  errors  were  from  H-optimal  followed  by  S^-optimal  for  (pi,P2,P3),  and  for  p±  it  was  ID- 
optimal  followed  by  S'H-optimal.  For  constraint  implementations  (C3)  and  (C'4)  (Tables  18  and  19), 
SlE-optimal  followed  by  H-optimal  had  the  smallest  standard  errors  for  (pi,p2,p:$),  and  ID-optimal 
followed  by  H-optimal  had  the  smallest  standard  errors  for  p^. 

In  conclusion,  ID-optimal  tended  to  have  the  smallest  standard  errors  for  p4,  whereas  S'F’-optimal 
or  H-optimal  had  the  smallest  standard  errors  for  (p \,p2, P‘3 ) ■  In  the  next  section  we  compute  the 
estimated  standard  errors  from  simulated  data  using  asymptotic  theory  and  bootstrapping  as  a 
different  method  of  comparing  the  optimal  design  methods. 
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Table  19:  Approximate  asymptotic  standard  errors  from  the  asymptotic  theory  (14)  resulting  from  different  optimal 
design  methods  for  do  =  (pi,P2,P3,P4),  optimization,  using  SolvOpt,  with  constraint  implementation  (C 4). 


Method 

SE{p  i) 

SE(p2 ) 

SE(p3) 

SE(p4) 

S'C-optimal 

4.921  x  10~3 

6.995  x  10~3 

3.633  x  10-6 

4.796  x  10-4 

H-optimal 

8.767  x  10~3 

1.249  x  10~2 

6.405  x  10"6 

1.965  x  10"4 

C-optimal 

7.154  x  10“3 

1.020  x  10~2 

5.253  x  10-6 

2.302  x  10-4 

6.5  Result  for  the  Glucose  Regulation  Model  with  the  Inverse  Problem 

As  in  the  harmonic  oscillator  example,  we  use  the  inverse  problem  with  the  OLS  formulation 
to  obtain  parameter  estimates  and  standard  errors  from  both  asymptotic  theory  (15)  and  the  boot¬ 
strapping  method  (18).  We  create  simulated  noisy  data  corresponding  to  the  optimal  time  meshes 
(presented  in  the  previous  section)  in  agreement  with  our  statistical  model  (absolute  error,  with 
independent  error  processes  for  G  and  I)  assuming  true  values  6q  to  be  the  parameter  values  found 
in  Table  15  and  iid  noise  with  £j  ~  A/"(0,  Ug).  We  assume  the  true  variances:  cig  G  =  25  and  cxg  7  =  4. 
In  this  section  we  only  estimate  a  subset  of  the  parameters  6  =  {pipP2-,P3iPa)-  hi  addition  to  the 
estimates  and  standard  errors,  we  also  report  the  estimated  covariance  between  estimated  parame¬ 
ters  according  to  asymptotic  theory  (15)  and  bootstrapping  (18).  For  comparison  purposes  we  also 
present  these  results  for  a  uniform  grid  using  the  same  T  =  150  and  N  =  30. 

The  optimal  time  points  for  each  of  the  three  optimal  design  methods  are  the  same  as  computed 
in  the  previous  results  section,  and  are  plotted  with  the  model  in  Figs.  13-16  for  the  four  differ¬ 
ent  constraints.  The  parameter  estimates,  standard  errors  and  covariances  are  estimated  from  the 
asymptotic  theory  (15)  corresponding  to  these  optimal  meshes  are  given  in  Tables  20,  22,  24,  and 
26,  respectively  for  the  four  different  constraints.  The  parameter  estimates,  standard  errors,  and 
covariance  between  parameters  are  estimated  from  the  bootstrapping  method  (18)  corresponding 
to  these  optimal  meshes  are  given  in  Tables  21,  23,  25,  and  27,  respectively  for  the  four  different 
constraints.  In  each  of  the  tables  are  also  results  on  the  uniform  grid  of  time  points. 

6.6  Discussion  for  the  Glucose  Regulation  Model  with  the  Inverse  Problem 

Comparing  the  resulting  parameter  estimates  from  simulated  data  on  the  different  optimal 
meshes  to  the  true  parameter  values,  $g  =  (pi,P2,P3,P4)  =  (2.6  x  10”2, 2.5  x  10~2, 1.25  x  10-5,4.1  x 
1CP3),  we  find  there  is  no  optimal  design  method  that  is  always  favorable.  Using  either  asymptotic 
theory  or  bootstrapping  to  compute  parameter  estimates  for  different  optimal  design  methods  and 
different  constraints,  we  examine  how  close  the  parameter  estimates  are  to  the  true  values.  Often 
(but  not  always)  these  parameter  estimates  from  the  different  optimal  meshes  are  the  same  order  of 
magnitude  as  the  true  values. 

The  results  for  the  uniform  mesh  are  given  for  comparison.  In  most  cases,  the  optimal  design 
methods  produce  closer  parameter  estimates  with  smaller  standard  errors  and  covariances  (as  esti¬ 
mated  by  asymptotic  theory  and  bootstrapping)  than  the  uniform  mesh. 

Asymptotic  theory:  parameter  estimates. 

For  the  constraint  implementation  (Cl)  using  asymptotic  theory  (Table  20),  the  estimates  for 
Pi-,  P2,  P3,  and  p4  are  closest  to  the  true  values  for  D-optimal  followed  by  C-optimal.  In  other 
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Table  20:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  asymptotic  theory  (15)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,P2,P3,P4)  =  (2.6  x  1  CP2,  2.5  x 
10  2 ,  1.25  x  10_5,4.1  x  1CP3)  and  N  =  30,  optimization,  using  fmincon,  with  constraint  implementation  (Cl). 


SID-optimal 

ID-optimal 

ID-optimal 

Uniform 

Pi 

3.036  x  10~2 

2.303  x  10-2 

2.267  x  10~2 

2.045  x  10-2 

SE(pi) 

4.977  x  10"3 

1.138  x  10"2 

4.634  x  10~3 

1.056  x  10"2 

P2 

2.180  x  10“2 

2.723  x  10"2 

2.818  x  10~2 

3.607  x  10"2 

SE(p2) 

7.657  x  10”3 

1.660  x  10-2 

6.818  x  10-3 

1.536  x  10-2 

P3 

9.213  x  10”6 

1.414  x  10"5 

1.565  x  10~5 

1.766  x  10“5 

SE(P3) 

3.946  x  10"6 

8.421  x  10-6 

3.732  x  10-6 

7.787  x  10-6 

Pa 

3.544  x  10“3 

4.174  x  10“3 

4.238  x  10~3 

4.027  x  10“3 

SE(pa) 

7.822  x  10"4 

4.510  x  10"4 

1.140  x  10~3 

4.817  x  10"4 

Cov(pi,p2) 

-3.377  x  10"5 

-1.846  x  10-4 

-2.779  x  10~5 

-1.579  x  10“4 

Cov(pi,p3) 

-1.873  x  10"8 

-9.520  x  10“8 

-1.630  x  10"8 

-8.160  x  10“8 

Cov(pi,j54) 

5.458  x  10“7 

8.794  x  10-7 

1.117  x  10~6 

8.615  x  10-7 

Cov(p2,h) 

2.815  x  10”8 

1.383  x  10"7 

2.289  x  10~8 

1.181  x  10"7 

Cov(p2,P4) 

2.851  x  10”7 

-6.379  x  10-7 

2.679  x  10~7 

-5.341  x  10~7 

Cov(j53,  j34) 

-5.0551  x  10-10 

-5.785  x  10-10 

-1.308  x  10~9 

-5.605  x  10_1° 

constraint  implementations,  which  optimal  sampling  distribution  produced  estimates  closest  to  the 
true  values  was  different  depending  on  the  parameter. 

For  constraint  implementation  (C 2)  (Table  22),  parameters  estimates  of  ( pi,P2,P3 )  were  closest 
to  the  true  values  for  the  optimal  sampling  distributions  from  ID-optimal  and  then  SID-optimal.  For 
P4  the  closest  parameter  estimates  were  from  the  uniform  mesh,  followed  by  ID-optimal. 

For  constraint  implementation  (C 3)  (Table  24),  the  closest  parameter  estimate  for  p\  came 
from  ID-optimal  followed  by  SID-optimal.  For  {p2,P3,Pa)  the  closest  estimates  came  from  ID-optimal 
followed  by  STD-optimal  (for  P2,P3)  and  ID-optimal  (for  p4). 

For  the  last  constraint  implementation  (C 4)  (Table  26),  ID-optimal  followed  by  SID-optimal  had 
parameter  estimates  closest  to  the  true  values  for  parameters  (p\,p2,p:i)-  ID-optimal  followed  by 
SID-optimal  had  the  closest  estimate  of  p4. 

Asymptotic  theory:  standard  errors. 

Here  we  compare  the  optimal  design  methods  based  on  which  has  the  smallest  standard  error 
estimates.  Again,  the  results  are  dependent  on  the  parameter  and  the  constraint  implementation. 

For  the  first  constraint  implementation  (Cl)  (Table  20),  the  smallest  standard  errors  estimated 
using  asymptotic  theory  are  from  ID-optimal  followed  by  SID-optimal  for  parameters  (pi,P2,P3),  and 
followed  by  ID-optimal  for  p4. 

For  the  constraint  implementation  (C2)  (Table  22),  the  smallest  standard  error  for  parame¬ 
ters  (pi,P2,P3)  come  from  SID-optimal  followed  by  ID-optimal.  For  p4,  the  smallest  standard  error 
estimates  are  from  the  uniform  mesh  followed  by  ID-optimal. 

For  the  constraint  implementation  (C 3)  (Table  24),  the  standard  error  estimates  for  parameters 
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Table  21:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  bootstrap  method  (18)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,P2,P3,P4)  =  (2.6  x  1  CP2, 2.5  x 
10  2 ,  1.25  x  1CP5,4.1  x  1CP3),  M  =  1000  bootstraps  and  N  =  30,  optimization,  using  fmincon,  with  constraint 
implementation  (Cl). 


SD-optimal 

D-optimal 

D-optimal 

Uniform 

Pi 

2.908  x  10"2 

2.220  x  10"2 

2.073  x  10"2 

1.973  x  10~2 

SE(Pi) 

6.215  x  10"3 

8.052  x  10-3 

5.708  x  10“3 

8.563  x  10"3 

P2 

2.486  x  10"2 

2.855  x  10"2 

3.126  x  10"2 

3.730  x  10~2 

SE(P2) 

1.179  x  10-2 

1.169  x  10-2 

9.810  x  10-3 

1.279  x  10-2 

P3 

1.075  x  10"5 

1.552  x  10-5 

1.864  x  10-5 

1.916  x  10-5 

SE(P3) 

6.541  x  10"6 

6.570  x  10"6 

6.302  x  10"6 

8.017  x  10~6 

Pi 

3.688  x  10"3 

4.215  x  10-3 

3.809  x  10“3 

3.984  x  10-3 

SE(p4 ) 

3.743  x  10~4 

1.855  x  10"4 

6.223  x  10"4 

2.098  x  10”4 

Cov(pi,p2) 

-6.799  x  10~5 

-9.116  x  10~5 

-5.323  x  10-5 

-1.053  x  10~4 

Cov(j3x ,  /33) 

-3.868  x  10”8 

-5.198  x  10~8 

-3.479  x  10-8 

-6.722  x  10”8 

Cov(pi,p4) 

-3.337  x  10”7 

2.268  x  10"7 

1.075  x  10"7 

5.716  x  10~8 

Cov(p2,P3) 

7.452  x  10~8 

7.529  x  10-8 

6.005  x  10-8 

9.990  x  10~8 

Cov(p2 )  Pi) 

1.050  x  10"6 

-1.262  x  10"7 

7.158  x  10-7 

2.310  x  10-7 

Cov(p3,p4) 

5.432  x  10-10 

-8.735  x  10-11 

-1.465  x  10-10 

8.308  x  10_n 

(pi,P2,P3)  are  smallest  using  the  mesh  from  SD-optimal,  followed  by  ID-optimal  (for  p\ )  and  D- 
optimal  (for  P2,P3)-  For  parameter  p 4,  the  smallest  standard  error  is  from  D-optimal  followed  by 
D-optimal. 

For  the  last  constraint  implementation  (C4)  (Table  26),  D-optimal  has  the  smallest  standard 
errors  for  parameters  {pi,P31Pi),  followed  by  SD-optimal  (for  pi,P3)  and  ID-optimal  (for  p 4).  For  P2, 
the  smallest  standard  errors  are  from  SD-optimal  followed  by  D-optimal. 

Asymptotic  theory:  covariance  estimates. 

We  also  compare  the  optimal  design  methods  based  on  which  has  the  smallest  covariance  esti¬ 
mates  in  absolute  value. 

For  the  first  constraint  implementation  (Ci)  (Table  20),  the  smallest  in  absolute  value  covariance 
estimates  come  from  either  S'D-optimal  or  D-optimal  for  different  pairs  of  parameters. 

For  constraint  implementation  ( C2 )  (Table  22),  SD-optimal  or  D-optimal  have  the  smallest  in 
absolute  value  covariance  estimates. 

For  constraint  implementation  (C3)  (Table  24),  SD-optimal  or  D-optimal  have  the  smallest  in 
absolute  value  covariance  estimates,  except  for  Cov(p2,P4)  where  D-optimal  is  the  smallest. 

For  constraint  implementation  (C4)  (Table  26),  D-optimal  has  the  smallest  in  absolute  value 
covariance  estimates,  except  for  Cov(]5i,p2)  where  S'D-optimal  is  the  smallest. 

Bootstrapping:  parameter  estimates. 

Here  we  compare  the  optimal  design  methods  based  on  which  had  bootstrapping  parameter 
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Table  22:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  asymptotic  theory  (15)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,P2,P3,P4)  =  (2.6  x  1  CP2,  2.5  x 
10  2 ,  1.25  x  10_5,4.1  x  1CP3)  and  N  =  30,  optimization,  using  fmincon,  with  constraint  implementation  (C2). 


Sid-optinral 

id-optimal 

id-optimal 

Uniform 

Pi 

2.118  x  10"2 

2.232  x  10-2 

2.116  x  10~2 

2.045  x  10-2 

SE(pi) 

5.063  x  10"3 

8.596  x  10“3 

5.298  x  10~3 

1.056  x  10“2 

P2 

3.509  x  10"2 

3.337  x  HT2 

4.356  x  10~2 

3.607  x  10“2 

SE{P2) 

8.020  x  10"3 

1.139  x  10-2 

8.465  x  10-3 

1.536  x  10-2 

P3 

1.772  x  10"5 

1.628  x  10“5 

1.958  x  10~5 

1.766  x  10“5 

SE{P3) 

4.247  x  10"6 

6.573  x  10-6 

4.874  x  10-6 

7.787  x  lO"6 

P4 

4.486  x  10"3 

3.993  x  10“3 

4.249  x  10~3 

4.027  x  10-3 

SE{pa) 

9.537  x  10"4 

5.919  x  10"4 

1.607  x  10~3 

4.817  x  10“4 

Cov(pi,p2) 

-3.569  x  10“5 

-9.416  x  10"5 

-3.811  x  10-5 

-1.579  x  10“4 

Cov(/)i  ,7)3) 

-2.036  x  10“8 

-5.566  x  10“8 

-2.376  x  10“8 

-8.160  x  10“8 

Cov(pi,p4) 

6.620  x  10"7 

1.227  x  10"7 

1.774  x  10-6 

8.615  x  10“7 

Cov(p2,P3) 

3.131  x  10"8 

7.280  x  10“8 

3.585  x  10~8 

1.181  x  10-7 

Cov(p2,p4) 

4.626  x  10"7 

4.238  x  10“7 

9.670  x  10~7 

-5.341  x  lO"7 

Cov(p3,7)4) 

-6.824  x  10~10 

9.532  x  10-13 

-2.353  x  10“9 

-5.605  x  10~10 

estimates  closest  to  the  true  values.  Often  these  results  are  different  for  the  different  parameters,  as 
well  as  the  constraint  implementation. 

For  the  first  constraint  implementation  (C\)  (Table  21),  bootstrapping  parameter  estimates  for 
(p\,p2-/P:s)  were  closest  to  the  true  values  for  Sid-optimal  followed  by  id-optimal.  For  p4,  the  closest 
parameter  estimates  came  from  id-optimal  and  then  FI-optimal. 

For  constraint  implementation  (C'2)  (Table  23),  parameter  estimates  for  p\  the  closest  parameter 
estimates  came  from  id-optimal  followed  by  Tl-optimal.  For  p2,  the  closest  parameter  estimates  came 
from  the  uniform  mesh  followed  by  id-optimal.  For  p 3,  the  uniform  mesh  then  S' FI-optimal  had  the 
closest  parameter  estimates  to  the  true  value.  For  p 4  the  closest  estimate  came  from  id-optimal 
followed  by  id-optimal. 

For  constraint  implementation  (C 3)  (Table  25),  the  closest  estimates  for  p\  came  from  id-optimal 
followed  by  id-optimal.  For  ( P2,P3,P4 )  the  closest  estimates  came  from  id-optimal  followed  by  SE- 
optimal  (for  752,7*3)  and  id-optimal  (for  754). 

For  the  last  constraint  (C4)  (Table  27),  the  closest  estimates  for  {pi,p2,P3)  came  from  id-optimal 
and  then  Sid-optimal.  For  p4  the  closest  estimate  to  the  true  value  came  from  Sid-optimal  followed 
by  id-optimal. 

None  of  the  optimal  design  methods  are  consistent  with  parameter  estimates  that  are  the  closest 
to  the  true  values  for  all  cases. 

Bootstrapping:  standard  errors. 

We  compare  the  optimal  design  methods  based  on  how  small  their  standard  errors  are  as  esti¬ 
mated  by  the  bootstrap  method. 
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Table  23:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  bootstrap  method  (18)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,P2,P3,P4)  =  (2.6  x  1  CP2, 2.5  x 
10  2 ,  1.25  x  1CP5,4.1  x  1CP3),  M  =  1000  bootstraps  and  N  =  30,  optimization,  using  fmincon,  with  constraint 
implementation  (C 2). 


SC-optimal 

C-optimal 

C-optimal 

Uniform 

Pi 

1.874  x  10~2 

1.883  x  10"2 

2.104  x  10"2 

1.973  x  10"2 

SE(pi) 

6.619  x  10-3 

8.291  x  10-3 

6.397  x  10“3 

8.563  x  10“3 

P2 

4.034  x  10~2 

4.249  x  10"2 

4.337  x  10"2 

3.730  x  10"2 

SE(P2 ) 

1.305  x  10-2 

1.458  x  10-2 

1.409  x  10-2 

1.279  x  10-2 

fa 

2.124  x  10~5 

2.069  x  10-5 

2.075  x  10-5 

1.916  x  10-5 

SE(P  s) 

8.241  x  10~6 

8.799  x  10"6 

7.733  x  10“6 

8.017  x  10"6 

fa 

4.341  x  10~3 

3.920  x  10-3 

3.988  x  10“3 

3.984  x  10-3 

SE(fa) 

4.228  x  10~4 

3.107  x  10-4 

6.192  x  10"4 

2.098  x  10"4 

Cov(pi,p2) 

-8.157  x  10-5 

-1.149  x  10~4 

-7.952  x  10-5 

-1.053  x  10~4 

Cov(/3 1,/)3) 

-5.272  x  10-8 

-7.128  x  10~8 

-4.687  x  10“8 

-6.722  x  10“8 

Cov(pi,p4) 

1.240  x  10"8 

1.275  x  10"7 

-5.657  x  10"7 

5.716  x  10“8 

Cov(p2  j  P3) 

1.048  x  10~7 

1.249  x  10-7 

1.042  x  Hr7 

9.990  x  10-8 

Cov(j)2,p4) 

4.220  x  10"7 

8.311  x  10“8 

2.226  x  10"6 

2.310  x  10"7 

Cov(p3,p4) 

6.133  x  10~n 

-2.764  x  10”11 

9.390  x  10~10 

8.308  x  10-n 

Comparing  the  standard  error  estimates  from  the  first  constraint  implementation  (Cl)  (Table  21) 
we  find  that  for  parameters  (pi,P2, P:i )  the  optimal  mesh  from  C-optimal  has  the  smallest  standard 
errors  followed  by  SC-optimal.  For  p4,  the  smallest  standard  errors  come  from  C-optimal  followed 
by  C-optimal. 

For  the  second  constraint  implementation  (C2)  (Table  23),  the  smallest  standard  errors  for 
parameters  (p\ ,  p;s)  are  from  C-optimal  followed  by  SC-optimal.  For  p2 ,  the  uniform  mesh  has 
the  smallest  standard  errors,  followed  by  SC-optimal.  For  p4,  the  uniform  mesh  has  the  smallest 
standard  errors  followed  by  the  optimal  mesh  from  C-optimal. 

For  the  constraint  implementation  (C3)  (Table  25),  the  smallest  standard  errors  come  from 
SC-optimal  for  parameters  (p\  ,p2)  followed  by  C-optimal.  For  parameters  (p3,p4)  the  optimal  mesh 
from  C-optimal  has  the  smallest  standard  errors,  followed  by  SC-optimal  (for  p%)  and  C-optimal 
(for  p4). 

For  the  last  constraint  implementation  (C4)  (Table  27),  the  smallest  standard  errors  for  pa¬ 
rameters  (pi,P2,P3)  are  from  C-optimal  followed  by  SC-optimal.  For  parameter  p4,  the  smallest 
standard  errors  are  from  C-optimal  followed  by  C-optimal. 

Bootstrapping:  covariance  estimates. 

For  the  first  constraint  implementation  (Cl)  (Table  21),  the  smallest  in  absolute  value  covariance 
estimates  as  estimated  by  the  bootstrapping  method  came  from  the  optimal  meshes  of  C-optimal 
(for  Cov(p2,Pi))  or  C-optimal  (for  Cov(])i,p2),  Cov(pi,p3),  and  Cov(p2lP3))  or  the  uniform  mesh 
(for  Cov(pi,p4)  and  Co v(p3,p4)). 
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Table  24:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  asymptotic  theory  (15)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,P2,P3,P4)  =  (2.6  x  1  CP2,  2.5  x 
10  2 ,  1.25  x  10_5,4.1  x  1CP3)  and  N  =  30,  optimization,  using  fmincon,  with  constraint  implementation  (C3). 


S'F'-optimal 

H-optimal 

Fi-optimal 

Uniform 

Pi 

2.960  x  10"2 

2.851  x  10-2 

2.970  x  10"2 

2.045  x  10"2 

SE(Pi) 

5.210  x  10"3 

8.937  x  10“3 

9.018  x  10"3 

1.056  x  10"2 

P2 

1.894  x  10"2 

1.303  x  10"2 

1.951  x  10"2 

3.607  x  10"2 

SE(P2) 

7.987  x  10"3 

1.853  x  10"2 

1.338  x  10"2 

1.536  x  10"2 

P3 

9.558  x  10"6 

8.981  x  10"6 

1.018  x  10"5 

1.766  x  10"5 

SE(h) 

4.201  x  10-6 

7.459  x  10-6 

6.701  x  10"6 

7.787  x  10"6 

Pa 

3.915  x  10"3 

3.945  x  10“3 

4.166  x  10"3 

4.027  x  10"3 

SE(Pa) 

8.373  x  10“4 

2.882  x  10"4 

4.366  x  10"4 

4.817  x  10~4 

Cov(pi,j52) 

-3.688  x  10~5 

-1.563  x  10~4 

-1.168  x  10-4 

-1.579  x  10-4 

Cov(pi,j53) 

-2.081  x  10~8 

-6.595  x  10"8 

-5.988  x  10“8 

-8.160  x  10~8 

Cov(pi,p4) 

5.971  x  10-7 

1.045  x  10-8 

5.570  x  10-7 

8.615  x  10-7 

Cov(p2,P3) 

3.113  x  10"8 

1.353  x  HT7 

8.834  x  10"8 

1.181  x  10"7 

Cov(p2,P4) 

3.520  x  10"7 

4.275  x  10-7 

-2.100  x  10-7 

-5.341  x  10~7 

Cov(p3,p4) 

-5.579  x  10'10 

5.546  x  10-n 

-3.461  x  10~10 

-5.605  x  10”10 

For  constraint  implementation  (C 2)  (Table  23),  D-optimal  or  U-optimal  have  the  smallest  in 
absolute  value  covariance  estimates,  except  for  Cov(pi,p4)  where  5'T-optimal  is  the  smallest  and 
Cov(p2 ,  P'i)  where  the  uniform  mesh  is  the  smallest. 

For  the  third  constraint  implementation  {C 3)  (Table  25),  Zl-optimal  or  S'£,-optimal  have  the 
smallest  in  absolute  value  covariance  estimates,  except  for  Cov{p21Pa)  where  U-optimal  is  the  small¬ 
est. 

For  the  last  constraint  implementation  ((74)  (Table  27),  the  smallest  in  absolute  value  covariance 
estimates  are  from  Fi-optimal,  except  for  Cov(pi,p4)  where  SU-optimal  is  the  smallest. 

Comparing  the  optimal  design  methods  based  on  the  bootstrapping  covariance  estimates,  we 
find  there  is  not  one  method  that  is  always  favorable. 
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Table  25:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  bootstrap  method  (18)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,P2,P3,P4)  =  (2.6  x  IIP2, 2.5  x 
10  2 ,  1.25  x  1(P5,4.1  x  10— 3 ) ,  M  =  1000  bootstraps  and  N  =  30,  optimization,  using  fmincon,  with  constraint 
implementation  (C3). 


S'U-optimal 

H-optimal 

i?-optimal 

Uniform 

Pi 

2.812  x  10~2 

2.738  x  10"2 

2.749  x  10~2 

1.973  x  10"2 

5.495  x  10~3 

6.344  x  10“3 

7.838  x  10~3 

8.563  x  10“3 

f>2 

2.160  x  10~2 

1.513  x  10-2 

2.300  x  10-2 

3.730  x  10-2 

9.465  x  10~3 

1.166  x  10"2 

1.183  x  10~2 

1.279  x  10"2 

P3 

1.122  x  10"5 

1.030  x  10-5 

1.240  x  10-5 

1.916  x  10-5 

5.360  x  10~6 

5.291  x  10"6 

6.429  x  10~6 

8.017  x  lO"6 

P4 

3.695  x  10"3 

4.011  x  10“3 

4.188  x  10~3 

3.984  x  10“3 

SE(pa) 

3.244  x  10“4 

1.311  x  10"4 

1.946  x  10~4 

2.098  x  10-4 

-4.787  x  10"5 

-6.858  x  10”5 

-8.971  x  10“5 

-1.053  x  10“4 

Cov(pi,p3) 

-2.768  x  10“8 

-3.288  x  10”8 

-4.933  x  10-8 

-6.722  x  10"8 

-1.523  x  10"7 

2.531  x  HT8 

9.876  x  10~8 

5.716  x  10“8 

Cov(p2,p3) 

4.913  x  10~8 

5.882  x  10-8 

7.443  x  10-8 

9.990  x  10-8 

Cov(p2  j  Pa) 

6.607  x  10~7 

9.576  x  10-8 

9.573  x  10-8 

2.310  x  10-7 

3.331  x  10”10 

3.083  x  10"11 

4.247  x  10”11 

8.308  x  10"11 
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Table  26:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  asymptotic  theory  (15)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,p2,P3,P4)  =  (2.6  x  1  CP2,  2.5  x 
10  2 ,  1.25  x  10_5,4.1  x  1CP3)  and  N  =  30,  optimization,  using  fmincon,  with  constraint  implementation  (C4). 


S'E'-optimal 

H-optimal 

i?-optimal 

Uniform 

Pi 

3.285  x  10"2 

2.540  x  10"2 

3.401  x  10"2 

2.045  x  10"2 

SE(Pi) 

6.396  x  10-3 

1.102  x  10-2 

6.353  x  10"3 

1.056  x  10"2 

P2 

1.783  x  10"2 

2.452  x  10"2 

1.079  x  10"2 

3.607  x  10"2 

SE(P2) 

8.983  x  10“3 

1.562  x  10"2 

9.052  x  10"3 

1.536  x  10"2 

P3 

9.094  x  10-6 

1.226  x  10-5 

6.025  x  10-6 

1.766  x  10"5 

SE(h ) 

5.375  x  10"6 

8.277  x  10“6 

4.850  x  10"6 

7.787  x  10"6 

Pa 

3.980  x  10-3 

3.958  x  10-3 

4.040  x  10"3 

4.027  x  10"3 

SE(Pi ) 

1.362  x  10“3 

4.896  x  10"4 

4.257  x  10~4 

4.817  x  10~4 

Cov(pi,p2) 

-5.251  x  10”5 

-1.692  x  10"4 

-5.587  x  10“5 

-1.579  x  10~4 

Cov(pi,ps) 

-3.260  x  10”8 

-9.060  x  10“8 

-3.043  x  10-8 

-8.160  x  10~8 

Cov(p4,p4) 

2.000  x  10"6 

1.1481  x  10"6 

4.352  x  10"7 

8.615  x  10"7 

Cov{p2,h) 

4.436  x  10-8 

1.281  x  10-7 

4.316  x  10"8 

1.181  x  10"7 

Cov(p2,P4) 

-2.356  x  10”7 

-1.016  x  10"6 

-2.022  x  10"7 

-5.341  x  10~7 

Cov(p3,p4) 

-2.203  x  10"9 

-8.152  x  10"10 

-3.102  x  10~10 

-5.605  x  10”10 
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Table  27:  Estimates,  standard  errors,  and  covariances  between  parameters  from  the  bootstrap  method  (18)  resulting 
from  different  optimal  design  methods  (as  well  as  for  the  uniform  mesh)  for  do  =  (pi,P2,P3,P4)  =  (2.6  x  IIP2, 2.5  x 
10  2 ,  1.25  x  IIP5, 4.1  x  10— 3 ) ,  M  =  1000  bootstraps  and  N  =  30,  optimization,  using  fmincon,  with  constraint 
implementation  (C4). 


S'E'-optimal 

H-optimal 

i?-optimal 

Uniform 

Pi 

3.243  x  10"2 

2.409  x  10"2 

3.304  x  10"2 

1.973  x  HP2 

4.851  x  10~3 

7.282  x  10“3 

3.410  x  10"3 

8.563  x  10“3 

f>2 

1.899  x  10~2 

2.659  x  10-2 

1.262  x  10"2 

3.730  x  10"2 

6.864  x  10~3 

1.056  x  10"2 

5.304  x  10"3 

1.279  x  HP2 

P3 

9.688  x  10~6 

1.385  x  10-5 

6.571  x  10-6 

1.916  x  10-5 

3.843  x  10~6 

5.780  x  10"6 

2.220  x  10"6 

8.017  x  10"6 

PA 

4.117  x  10"3 

3.967  x  10“3 

4.025  x  10"3 

3.984  x  10“3 

SE(pa ) 

7.050  x  10~4 

2.147  x  10"4 

1.892  x  10-4 

2.098  x  10-4 

-2.970  x  10“5 

-7.462  x  10”5 

-1.700  x  10“5 

-1.053  x  10~4 

Cov(pi,p3) 

-1.680  x  10-8 

-4.139  x  10~8 

-7.230  x  10-9 

-6.722  x  10”8 

2.879  x  10~8 

6.162  x  10“8 

9.816  x  10"8 

5.716  x  HP8 

Cov(p2,P3 ) 

2.477  x  10~8 

5.986  x  10-8 

1.133  x  10"8 

9.990  x  10-8 

Cov(p2,Pa) 

6.321  x  10-7 

1.286  x  10-7 

2.627  x  10"9 

2.310  x  HP7 

4.890  x  10"11 

2.363  x  10"11 

-1.770  x  10"11 

8.308  x  10"11 
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7  Conclusions 


We  compared  ID-optimal,  H-optimal  and  .SH-optimal  design  methods  for  a  simple  differential 
equation  model:  the  logistic  population  model,  a  second  order  differential  equation:  the  harmonic 
oscillator  model,  and  a  vector  system  for  glucose  regulation.  ID-optimal  and  H-optimal  design  meth¬ 
ods  are  more  established  in  the  literature.  Our  comparisons  test  the  performance  of  SiH-optimal 
design,  which  is  a  relatively  newer  method. 

For  the  logistic  example,  the  optimal  design  methods  were  compared  using  the  Monte  Carlo 
method  for  asymptotic  standard  errors.  Comparing  the  average  and  median  parameter  estimates  to 
their  true  values,  we  find  that  .SH-optimal  has  closest  parameter  estimates  for  N  =  10  time  points. 
For  N  =  15,  no  method  had  estimates  that  were  always  closest  to  the  true  values.  In  all  cases 
each  optimal  design  methods  produced  estimates  close  to  the  true  values.  The  average  and  median 
standard  errors  for  I\  were  smallest  from  the  optimal  mesh  from  F-optimal.  For  parameters  r  and 
.To,  SlE-optimal  had  the  smallest  average  and  median  standard  errors.  Overall,  no  optimal  design 
method  is  consistently  favorable  for  this  logistic  example. 

For  the  harmonic  example,  comparing  the  approximate  asymptotic  standard  errors,  we  found 
that  different  optimal  design  methods  were  favorable  for  different  parameters.  ID-optimal  often  had 
the  smallest  standard  errors  for  K  and  x\.  AH-optimal  often  had  the  smallest  standard  errors  for 
C.  For  X2,  either  SFi-optimal  or  H-optimal  had  the  smallest  standard  errors.  We  also  compared 
methods  using  the  inverse  problem  with  simulated  data  and  asymptotic  theory  and  bootstrapping. 
Comparing  methods  based  on  who’s  parameter  estimates  were  closest  to  the  true  values,  and  who 
had  the  smallest  standard  errors  or  covariances,  there  was  no  method  that  was  preferred  over  the 
others.  In  each  comparison,  the  best  optimal  design  method  often  depended  on  the  constraint 
implementation,  the  choice  of  T  =  14.14  or  T  =  28.28,  and  the  parameter. 

For  the  glucose  regulation  model,  comparing  the  approximate  asymptotic  standard  errors,  we 
found  that  for  parameters  (p\ ,  p2-,P:i)  either  S'Fi-optimal  or  H-optimal  had  the  smallest  standard 
errors.  ID-optimal  tended  to  have  the  smallest  standard  errors  for  p±.  We  also  compared  the  optimal 
design  methods  for  the  inverse  problem  using  asymptotic  theory  and  bootstrapping.  Comparing 
the  parameter  estimates  to  their  true  values,  none  of  the  optimal  design  methods  were  consistently 
closer.  Comparing  the  optimal  design  methods  based  on  who  had  the  smallest  standard  errors 
and  covariances  we  found  that  no  method  was  preferable  over  the  others.  However,  the  optimal 
design  methods  often  had  smaller  standard  errors  and  covariances  than  the  uniform  mesh.  The 
constraint  implementation,  parameter,  and  choice  of  asymptotic  theory  or  bootstrapping  influenced 
which  optimal  design  method  would  be  favorable  for  this  example. 

The  best  choice  of  optimal  design  method  depends  on  the  complexity  of  the  model,  the  type  of 
constraint  one  is  using,  the  subset  of  parameters  you  are  estimating,  and  even  the  choice  of  N  and 
T.  The  examples  in  this  comparison  provide  evidence  that  S'F’-optimal  design  is  competitive  with 
ID-optimal  and  H-optimal  design,  and  in  some  cases  S'H-optimal  design  is  a  more  favorable  method. 
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Appendix-Constraints  and  Implementation  of  Constrained  Optimiza¬ 
tion 


We  used  several  constrained  optimization  algorithms  to  solve  the  grid  selection  minimization 
problem  of  the  form 

u*  =  ruin  J(u), 


subject  to  the  constraint (s) 


Au  <  b,  and/or  Aequ  =  beq, 


where  u  is  a  IV-vector,  A  is  a  (N  +  1)  x  N  matrix,  b  is  a  (TV  +  l)-vector,  Aeq  is  a  N  x  N  matrix,  and 
beq  is  a  scalar. 


For  our  problem,  we  have  the  constraint 

0  <  zzi  <  zz2  <  . . .  <  ZW  <  1, 
where  t  =  (t\,  -.-An)  =  vT  =  (u\T, ...,  unT ),  then 

0  <  ti  <t,2  <  ...  <tN  <T. 

To  express  this  constraint  in  the  form 

Au  <  b ,  and/or  Aequ  =  beq, 

we  have  several  options  in  algebraic  formulations.  Our  four  different  constraint  implementations  are 
detailed  below  and  the  differences  in  the  implementations  of  the  constrained  optimization  algorithm 
account  for  the  differences  in  the  optimal  meshes  generated.  As  is  explained,  a  primary  difference  in 
carrying  out  the  optimizations  is  the  number  of  points  over  which  we  optimize  (i.e.,  the  number  of 
degrees  of  freedom  in  the  problem). 


Constraint  implementation  (Cl):  For  this  constraint  implementation,  it  differs  from  the  other 
three  in  that  it  is  not  required  that  the  end  points  are  included  in  the  optimal  mesh.  For  this 
constraint  we  define  the  (N  +  1)  x  N  matrix, 


We  define  the  ( N  +  l)-vector 
The  constraint  Au  <  b,  implies 
Setting  t  =  uT,  we  obtain 


A  = 


/  1  0  0  0 
-110  0 
0-110 


\  0  •••  0 


-1  1 


b=  [0,  -  -  -  ,0,lf 


0  <  U\  <  U2  <  .  .  .  <  ZZ/v  —  1- 


o  <  ti  <  t2  <  ...  <  tN  <  T. 
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In  this  case  we  optimize  over  N  points. 

Constraint  implementation  ((72): 

For  this  constraint  implementation,  we  require  that  the  end  points  are  included  in  the  optimal 
mesh.  We  optimize  over  the  remaining  mesh  points  (f2,  •  • .  tjv- 1)-  For  this  constraint  we  define  the 
(IV  —  1)  x  ( N  —  2)  matrix, 


We  define  the  ( N  —  l)-vector 

b=  [0,  -  -  -  ,0,1]T. 

The  constraint  Av  <  b ,  implies 

0  =  V\  <  V-2  <  ^2  <  •  •  •  <  Z^AT-l  <  on  =  1  • 

Upon  setting  t  =  CT,  we  obtain 

0  =  t\  <  ^2  ^  •  •  •  5:  tjV-1  7  tjV  =  T. 

In  this  case  we  optimize  over  N  —  2  points. 

Constraint  implementation  ((73): 

For  the  third  constraint  implementation,  we  include  the  end  points  in  the  optimal  mesh.  We  optimize 
over  the  remaining  mesh  points  (£2,  •  •  •  tjv-i)-  For  this  constraint  we  define  the  (IV  —  1)  x  (IV  —  2) 
matrix, 


We  define  the  (IV  —  1  (-vector, 

b=[ 0,--  -  ,0,T]t. 

The  constraint  <  6,  implies 

u-j  >  0,  for  i  =  2, . . . ,  N  —  1 

and 

+  zaj  . . .  +  i/jv-i  <  T. 

To  form  t  from  v,  we  first  must  define  the  (IV  —  2)  x  (IV  —  2)  matrix 
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Setting  t\  =  0,  tw  =  T  and 


which  implies  that 


Then 


or  equivalently, 


[t2,  ■  ■  •  ,tjV-l]T  =  B[v2,  ,  vn-i]  r, 


tk  =  z/j,  for  all  k  =  2, . . .  TV  —  1. 


0  =  fl  <  V2  <  V2  +  ^3  <  •  •  •  <  (v2  +  ^3  +  •  •  •  +  VN- 1)  <  tv  =  T, 


0  =  ti  <  t2  <  •  •  •  <  tjv— i  £  t./v  =  T. 


We  again  optimize  over  TV  —  2  points. 

Constraint  implementation  ((74): 

For  the  fourth  constraint,  we  include  the  end  points  in  the  optimal  mesh.  For  this  constraint  we 
define  the  (A/"  —  1)  x  (N  —  1)  matrix, 

/  -1  0  0  0  •••  \ 

0-100  ••• 

A  =  0  0-10  ••• 


0  0-1 


We  define  the  (TV  —  l)-vector, 


The  constraint  Av  <  b,  implies 


b=\ 0,---  ,0]' 


Vi  >  0,  for  i  =  2, . . . ,  TV. 


In  addition,  we  define  the  (TV  —  l)-row  vector 

A-eq  =  [1,  1,  •  •  •  ,  1] , 

and  the  scalar  beq  =  T.  The  additional  constraint,  Aeqv  =  beq,  implies 


To  form  t  from  v.  we  first  must  define  the  (TV  —  1)  x  (TV  —  1)  matrix 

/  1  0  0  0  •••  \ 

11  0  0  ••• 

B=  1  1  1  0  •••  . 


1  •••  1  1  1 
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Setting  t\  =  0  and 


[t2,...,t  at]1  =  B[v2,.  .  .  ,VN]T, 


which  implies  that 


k 

tk  =  i 'j,  for  all  k  =  2, . . .  N. 

3= 2 


Then 


0  =  tl  <  ^2  <  ^2  +  ^3  <  •  •  •  <  (^2  +  ^3  +  •  •  •  +  Yv)  =tN  =  T, 


or  equivalently, 


0  =  ii  <  t2  <  . . .  <  _ l  £  i/v  =  T. 


In  this  algorithm  we  again  effectively  optimize  over  N  —  2  points. 
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